SetDirectory[$HomeDirectory<>"/Documents/handouts/math467/IterativeMethods"]

Gauss-Seidel

g[x_,0]:=Exp[x]

g[x_,1]:=Sin[π x/2]

g[0,y_]:=1-(y^2)

g[1,y_]:=Exp[1-y]

nsoln=First[NDSolve[{D[u[x,y],x,x]+D[u[x,y],y,y]==0,u[0,y]==g[0,y],u[x,1]==g[x,1],u[1,y]==g[1,y],u[x,0]==g[x,0]},u,{x,0,1},{y,0,1}]]

m=n=10;h=1/n

Table[-u[i+1,j]-u[i-1,j]+4u[i,j]-u[i,j+1]-u[i,j-1]==0,{j,1,m-1},{i,1,n-1}];

%/.{u[0,j_]->g[0,j h],u[i_,0]->g[i h,0],u[n,j_]->g[n h,j h],u[i_,m]->g[i h, m h]};

{b,a}=CoefficientArrays[Flatten[%],Flatten[Table[u[i,j],{j,1,m-1},{i,1,n-1}]]]

Normal[-b]//N

GSsoln=Inverse[a].(-b)//N

isoln=Join[{Table[g[i h,0],{i,0,n}]},Join[Transpose[{Table[g[0,j h],{j,1,m-1}]}],Partition[GSsoln,9,9],Transpose[{Table[g[n h,j h],{j,1,m-1}]}],2],{Table[g[i h,m h],{i,0,n}]}]//N

xs=Table[FindRoot[z==u[x,1-x]/.nsoln,{x,0.85}],{z,0.25,2.5,0.25}]

ContourPlot[u[x,y]/.nsoln,{x,0,1},{y,0,1},FrameLabel->{"x","y"},Contours->Range[0.25,2.5,0.25],Epilog->{Text[NumberForm[u[x,1-x]/.nsoln,{3,2}],{x,1-x},Background->White]/.xs}]

Export["exampleLaplace.pdf",%,"PDF"]

Export["exampleLaplace.eps",%%,"EPS"]

Successive Over-Relaxation (SOR) Method

a={{-4,1,0,0,0},{1,-4,1,0,0},{0,1,-4,1,0},{0,0,1,-4,1},{0,0,0,1,-4}};

MatrixForm[a]

Inverse[a].Range[5]

Inverse[a].Range[5]//N

dmat=DiagonalMatrix[Diagonal[a]];

lmat=LowerTriangularize[a,-1];

umat=UpperTriangularize[a,1];

tjmat=Inverse[a].(lmat+umat);

Eigenvalues[tjmat]

Max[Abs[Eigenvalues[tjmat]]]

2/(1+Sqrt[1-%*%])

N[%]

Example

Clear[x,y,nsoln,u]

nsoln=First[NDSolve[{D[u[x,y],x,x]+D[u[x,y],y,y]==Exp[x-y],u[0,y]==0,u[x,0]==0,u[x,3]==3x,u[4,y]==4y},u,{x,0,4},{y,0,3}]]

Clear[m,n,u,f, g,h]

h=1/10;m=30;n=40;

f[x_,y_]:=Exp[x-y]

g[x_,y_]:=x y

Table[20u[i,j]-4(u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1])-(u[i-1,j-1]+u[i-1,j+1]+u[i+1,j-1]+u[i+1,j+1])==-(h^2)(f[i-1,j]+f[i+1,j]+8f[i,j]+f[i,j-1]+f[i,j+1])/2,{j,1,m-1},{i,1,n-1}];

Table[20u[i,j]-4(u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1])-(u[i-1,j-1]+u[i-1,j+1]+u[i+1,j-1]+u[i+1,j+1])==-6(h^2)(f[i h,j h]+(h^2)(Derivative[2,0][f][i h,j h]+Derivative[0,2][f][i h,j h])/12),{j,1,m-1},{i,1,n-1}];

%/.{u[0,j_]->g[0,j h],u[n,j_]->g[n h,j h],u[i_,0]->g[i h,0],u[i_,m]->g[i h, m h]};

{b,a}=CoefficientArrays[Flatten[%],Flatten[Table[u[i,j],{j,1,m-1},{i,1,n-1}]]]

dmat=DiagonalMatrix[Diagonal[a]]

lmat=LowerTriangularize[a,-1]

umat=UpperTriangularize[a,1]

Max[Abs[Eigenvalues[N[Inverse[dmat].(lmat+umat)]]]]

2/(1+Sqrt[1-%*%])

DensityPlot[u[x,y]/.nsoln,{x,0,4},{y,0,3},FrameLabel->{"x","y"},PlotTheme->"Detailed",AspectRatio->Automatic]

xs=Table[FindRoot[z==u[x,3x/4]/.nsoln,{x,3.0}],{z,0.5,9.5,1.0}]

ContourPlot[u[x,y]/.nsoln,{x,0,4},{y,0,3},FrameLabel->{"\!\(\*TemplateBox[<|"boxes" -> FormBox[
StyleBox["x", "TI"], TraditionalForm], "errors" -> {}, "input" -> "x", "state" -> "Boxes"|>,"TeXAssistantTemplate"]\)","\!\(\*TemplateBox[<|"boxes" -> FormBox[
StyleBox["y", "TI"], TraditionalForm], "errors" -> {}, "input" -> "y", "state" -> "Boxes"|>,"TeXAssistantTemplate"]\)"},Contours->Range[0.5,9.5,1],AspectRatio->3/4,Epilog->{Text[u[x,3x/4]/.nsoln,{x,3x/4},Background->White]/.xs}]

Export["9ptPoissonSOR.pdf",%,"PDF"]

Export["9ptPoissonSOR.eps",%%,"EPS"]