Skip to main content Skip to footer site map
Reza Malek-Madani

parclap2html

(* This program solves the BVP laplace u = 0, with
boundary conditions u(x,0) = u(0, y) = 0, and 
u(x, 1) = f(x), u(1, y) = g(y).  To reduce the
time of execution when various f's and g's are entered,
it is best to compute the 4th order tensor A separately. *)

eps = 0.31;
enn=7; 
L[u_]:=D[u, {x, 2}] + D[u, {y, 2}];
newT[i_, x_] = ChebyshevT[i,2x -1]; 
T[k_, x_] =newT[k,x]+(newT[k,0]-newT[k,1])x - newT[k,0];
phi[k_,l_,x_,y_] = T[k,x] T[l, y];
m1[k_, l_, m_, n_]:=Integrate[L[phi[k,l,x,y]]*phi[m,n,x,y],
                 {x, 0, 1}, {y, 0, 1}];
A = Table[m1[k,l,m,n], {k,2,enn}, {l, 2, enn}, {m, 2, enn}, {n, 2, enn}];
f[x_] = x*(1-x); g[y_] = y*(1-y); u0[x,y]=f[x]*y+eps*g[y]*x;
ua[x_,y_] = u0[x,y]+Sum[a[m,n] phi[m,n,x,y], {m, 2, enn}, {n, 2, enn}];

m2[k_, l_] := Integrate[L[u0[x,y]] phi[k, l, x, y], {x,0,1},{y,0,1}];
B=Table[m2[k,l], {k,2,enn}, {l,2,enn}];
coefficients=Flatten[Table[a[m,n], {m, 2, enn}, {n, 2, enn}]];
equations=Flatten[Table[Sum[A[[m,n,k,l]] a[m+1,n+1], {m,1,enn-1}, 
      {n,1,enn-1}] == - B[[k,l]], {k,1,enn-1}, {l,1,enn-1}]];
sol=Solve[equations,coefficients];
u[x_,y_]=ua[x,y]/. sol;
graph1=ContourPlot[First[u[x,y]], {x,0,1},{y,0,1}, ColorFunction->Hue,
PlotLabel->"Galerkin with 7 Chebyshev Polynomials"];
(* Exact Solution  *)
a[n_]:= 2/Sinh[n Pi]*Integrate[f[x] Sin[n Pi x], {x, 0, 1}];
aa = Table[a[n], {n, 1, 20}];
exact1[x_, y_] = Sum[aa[[i]] Sin[i Pi x] Sinh[i Pi y], {i, 1, Length[aa]}];
b[n_]:= 2/Sinh[n Pi]*Integrate[eps*g[y] Sin[n Pi y], {y, 0, 1}];
bb = Table[b[n], {n, 1, 20}];
exact2[x_, y_] = Sum[bb[[i]] Sin[i Pi y] Sinh[i Pi x], {i, 1, Length[bb]}];
exact[x_,y_]=exact1[x,y]+exact2[x,y];
Print[N[exact[0.5,0.5]], " --- ", First[u[0.5,0.5]]]
graph2=ContourPlot[exact[x,y], {x,0,1},{y,0,1}, ColorFunction->Hue,
PlotLabel->"Contours of u -- 20 term Fourier series"];
(* Parcel Evolution *)
uxy=Simplify[Chop[First[u[x,y]]]];

Clear[output1,oldsolution,data,ColoredSnapshots1,diffeqn];
f[x_,y_] =0.1 D[uxy, y]; g[x_, y_] = -0.1 D[uxy, x];
m=1;                (*Number of parcels*)
dy = 0.73;           (*Distance between parcels at time 0*)
n=40;               (*Number of points on the boundary of parcels*)
diffeqn[tfinal_,a_,b_]:=NDSolve[{x'[t]==f[x[t],y[t]],
                                 y'[t]==g[x[t],y[t]],
                                 x[0]==a, y[0]==b},
                                 {x,y}, {t,0,tfinal}, MaxSteps->5000];
tfinal=60; snap=15;
oldsolution=Table[diffeqn[tfinal, 0.11+0.1*Cos[w], h*dy +0.1*Sin[w]],
                 {h, 1, m}, {w, 0, 2*Pi, 2*Pi/n}];
data=Table[{x[t], y[t]} /. oldsolution, {t, 0, tfinal, snap}];
ColoredSnapshots1=Table[Graphics[
              {RGBColor[0.2,0.3,0.5], 
         Map[Polygon, Flatten[data[[i]][[j]],1], {0}]}], {i, Length[data]},
         {j, 1, m}];
square=Graphics[{Line[{{0,0},{1,0}}], Line[{{1,0},{1,1}}],
                Line[{{1,1}, {0,1}}], Line[{{0,1}, {0,0}}]}];
output1=Show[square,ColoredSnapshots1, AspectRatio->Automatic,
PlotRange->All, PlotLabel->"Parcel motion using the Galerkin approximation"];
output=Show[graph1,output1];
go to Top