4620 Homework 6
Guide

Design Decisions

I need to solve the heat equation on a rectangular domain with peridiodic boundary conditions in the y direction and Neumann and Dirichlet in the x direction.  

Decisions:
1)  I am going to use Ghost Points and a symmetric difference for the Neumann Conditions. I am going to include the ghost points in the computation.
2) I am going to include the Dirichlet points in the computation.
3) I am going to implement the periodic conditions by including one set of the paired boundary conditions and I am going to use modular arithmetic to track the wrap around points.

Bookkeeping and Terminology

I am going to use the standard vertex labelling (or unrolling) that we have dicussed before.  The domain is 
	a≤x≤b
	c≤y≤d
I have nx and ny interior points in each direction. The coordinates and are labelled by the following functions

{{a,b},{c,d}}={{1.0,4.0},{2.0, 5.0}};
{nx, ny}={3,4};
{Δx,Δy}={((b-a)/(nx+1)),((d-c)/(ny+1))};
p[i_,j_]:={a+i Δx,c + j Δy}
Count1[i_,j_]:=1 + i + (nx+3)*j
ListPlot[
{
Flatten[Table[p[i,j],{i,1,nx+1},{j,1,ny}],1],
Table[p[0,j],{j,1,ny}],
Table[p[nx+2,j],{j,1,ny}],
Table[p[i,0],{i,1,nx+1}],
{p[0,0],p[nx+2,0]}
},
Prolog->{Pink,Rectangle[{a,c},{b,d}]},
Epilog->Table[ Text[Count1[i,j],p[i,j],{-1,-1}],
{i,0, nx+2},{j,0, ny}],
Frame->True,
Axes->False,
GridLines->{{a,b},{c,d}},
PlotRange->{{a-Δx,b+2 Δx},{c-Δy,d+Δy}},
PlotStyle->{Red, Blue, Green,Yellow,Cyan}
]

Periodic Counter

So I am going to make the counter work correctly for the top wrap around points.

Count2[i_,j_]:=1 + i + (nx+3)*Mod[j,ny+1]
ListPlot[
{
Flatten[Table[p[i,j],{i,1,nx+1},{j,1,ny}],1],
Table[p[0,j],{j,1,ny}],
Table[p[nx+2,j],{j,1,ny}],
Flatten[Table[p[i,j],{i,1,nx+1},{j,{0,ny+1}}],1],
{p[0,0],p[nx+2,0],p[0,ny+1],p[nx+2,ny+1]}
},
Prolog->{Pink,Rectangle[{a,c},{b,d}]},
Epilog->Table[ Text[Count2[i,j],p[i,j],{1,1}],
{i,0, nx+2},{j,0, ny+1}],
Frame->True,
Axes->False,
GridLines->{{a,b},{c,d}},
PlotRange->{{a-Δx,b+2 Δx},{c-Δy,d+Δy}},
PlotStyle->{Red, Blue, Green,Yellow, Cyan}
]

As long as I use Count2 I can treat all the yellow points as interior just like the red ones.  This is our periodic trick.

Equations

OK:
The red are interior points.
The blue are Dirichlet points.  
The green are ghost points. 
The yellow are periodic points that I can deal with using the modular trick.
The cyan points are peridic and ghost!

Equations: Interior and Periodic Points.

Each row in the matrix is an equation.  The equation for an interior point is 
	(Subscript[UNew, i,j])-(Subscript[UOld, i,j])=(1/2)((Δt)/((Δx)^2))((Subscript[UNew, i-1,j])-2(Subscript[UNew, i,j])+(Subscript[UNew, i+1,j]))+(1/2)((Δt)/((Δy)^2))((Subscript[UNew, i,j-1])-2(Subscript[UNew, i,j])+(Subscript[UNew, i,j+1]))+(1/2)((Δt)/((Δx)^2))((Subscript[UOld, i-1,j])-2(Subscript[UOld, i,j])+(Subscript[UOld, i+1,j]))+(1/2)((Δt)/((Δy)^2))((Subscript[UOld, i,j-1])-2(Subscript[UOld, i,j])+(Subscript[UOld, i,j+1])) + Δt (1/2)((Subscript[fOld, i,j])+(Subscript[fNew, i,j]))

Remember our notation UNew is the 2D array of values but uNew is the vector you get by squishing UNew flat.  Samething for UOld and uOld.  Moving UNew to one side and UOld to the other I get.
(1)	(Subscript[UNew, i,j])-(1/2)((Δt)/((Δx)^2))((Subscript[UNew, i-1,j])-2(Subscript[UNew, i,j])+(Subscript[UNew, i+1,j]))-(1/2)((Δt)/((Δy)^2))((Subscript[UNew, i,j-1])-2(Subscript[UNew, i,j])+(Subscript[UNew, i,j+1]))
on one side and all the stuff I know 
(2)	(Subscript[UOld, i,j])+(1/2)((Δt)/((Δx)^2))((Subscript[UOld, i-1,j])-2(Subscript[UOld, i,j])+(Subscript[UOld, i+1,j]))+(1/2)((Δt)/((Δy)^2))((Subscript[UOld, i,j-1])-2(Subscript[UOld, i,j])+(Subscript[UOld, i,j+1])) + ((Δt)/2)((Subscript[fOld, i,j])+(Subscript[fNew, i,j])) 
on the other.   I am going to build two matrices row by row.  Each matrix is square of size (ℝ^(n×n)) where for this problem n=(ny+1)(nx+3).

The matrix ANew is associated with the New side or  (1).  The matrix AOld is associated with the Old side or (2).  I can always keep track of where I am in the matrix using Count2[i,j].  If {i,j} then (1) is telling me the entires I need to insert into the Count2[i,j]th row of ANew and (2) is telling me what entries I need to insert into the same row of AOld.  I am going to call the rations {αx,αy}=(1/2){((Δt)/((Δx)^2)),((Δt)/((Δy)^2))}.  The following does exactly this in hopefully an organized and comprehensible manner.

InteriorPointMatrixBuild[{i_,j_,Count2_},{αx_,αy_},{AN_,AO_}] := Module[
{k=Count2[i,j],ANew=AN, AOld=AO},
ANew[[k,k]]=1+2αx + 2αy;
ANew[[k,Count2[i-1,j]]]=-αx;
ANew[[k,Count2[i+1,j]]]=-αx;
ANew[[k,Count2[i,j-1]]]=-αy;
ANew[[k,Count2[i,j+1]]]=-αy;
AOld[[k,k]]=1-2αx - 2αy;
AOld[[k,Count2[i-1,j]]]=αx;
AOld[[k,Count2[i+1,j]]]=αx;
AOld[[k,Count2[i,j-1]]]=αy;
AOld[[k,Count2[i,j+1]]]=αy;
{ANew,AOld}
]

Remember, this is the same equation at all the interior/periodic and because of the way we are doing it with Ghost points the "interior" Neumann points as  well.

Equations: Dirichlet Points.

Each row in the matrix is an equation.  The equation for a Dirichlet point is 
	(Subscript[UNew, i,j])= fD((Subscript[x, i]),(Subscript[y, j]),(Subscript[t, New])).
This is easy all we have to do is upgrade our ANew equation and do nothing to AOld.

The following does exactly this in hopefully an organized and comprehensible manner. It does not do anything with "f".

DirichletPointMatrixBuild[{i_,j_,Count2_},{αx_,αy_},{AN_,AO_}] := Module[
{k=Count2[i,j],ANew=AN, AOld=AO},
ANew[[k,k]]=1;
{ANew,AOld}
]

Equations: Neumann Ghost Points.

Each row in the matrix is an equation.  The equation for a Ghost point is 
	(Subscript[UNew, i,j])-(Subscript[UNew, i-2,j])=2Δx fN((Subscript[x, i-1]),(Subscript[y, j]),(Subscript[t, New])).
This is easy all we have to do is upgrade our ANew equation and do nothing to AOld.

The following does exactly this in hopefully an organized and comprehensible manner. It does not do anything with "fN".

GhostPointMatrixBuild[{i_,j_,Count2_},{αx_,αy_},{AN_,AO_}] := Module[
{k=Count2[i,j],ANew=AN, AOld=AO},
ANew[[k,k]]=1.0;
ANew[[k,Count2[i-2,j]]]=-1.0;
{ANew,AOld}
]

Forcing Terms.

We are setting up our iteration in the form
	ANew.uNew=AOld.uOld+(OverVector[f])
where (OverVector[f]) is defined in terms of our input fEq, fD, and fN depending on which sort of point we are at.   This should update the terms appropriately.

UpDateFInterior[{i_,j_,n_,Count2_},{Δx_,Δy_,Δt_,p_},{fEq_,fD_,fN_},f_]:= Module[
{k=Count2[i,j]},
f[[k]]=0.5*Δt(fEq[p[i,j],n*Δt]+fEq[p[i,j],(n+1)*Δt]);
f]
UpDateFDirichlet[{i_,j_,n_,Count2_},{Δx_,Δy_,Δt_,p_,Count2_},{fEq_,fD_,fN_},f_]:= Module[
{k=Count2[i,j]},
f[[k]]=fD[p[i,j],(n+1)*Δt];
f]
UpDateFNeumann[{i_,j_,n_,Count2_},{Δx_,Δy_,Δt_,p_,Count2_},{fEq_,fD_,fN_},f_]:= Module[
{k=Count2[i,j]},
f[[k]]=2*Δx*fN[p[i-1,j],(n+1)*Δt];
f]

Initial Conditions

We are setting up our iteration in the form
	ANew.uNew=AOld.uOld+(OverVector[f])
where (OverVector[f]) is defined in terms of our input fEq, fD, and fN depending on which sort of point we are at.   This should update the terms appropriately.

SetIC[{nx_, ny_,Count2_, p_},UO_]:= Module[
{U},
U=ConstantArray[0,(nx+3)(ny+1)];
Do[
U[[Count2[i,j]]]=UO[p[i,j]],
{i, 0, nx+2},
{j,0, ny}];
U
]

Putting it all together

First Lets build the matrices

{{a,b},{c,d}}={{1,4},{0.0, 2*N[Pi]}};
{nx, ny}={13, 14};
{Δx,Δy}={((b-a)/(nx+1)),((d-c)/(ny+1))};
p[i_,j_]:={a+i Δx,c + j Δy}
Count2[i_,j_]:=1 + i + (nx+3)*Mod[j,ny+1]
Δt=0.01;
n=(ny+1)(nx+3);
{αx,αy}=(1/2){((Δt)/((Δx)^2)),((Δt)/((Δy)^2))};
ANew=SparseArray[{},{n,n}];
AOld=SparseArray[{},{n,n}];

Do[
{ANew,AOld}=InteriorPointMatrixBuild[{i,j,Count2},{αx,αy},{ANew,AOld}],
{i,1,nx+1},
{j, 0, ny}
]
Do[
{ANew,AOld}=DirichletPointMatrixBuild[{0,j,Count2},{αx,αy},{ANew,AOld}],
{j, 0, ny}
]
Do[
{ANew,AOld}=GhostPointMatrixBuild[{nx+2,j,Count2},{αx,αy},{ANew,AOld}],
{j, 0, ny}
(* Here is a comment *)
]

{MatrixPlot[AOld],MatrixPlot[ANew]}

What deos the update scheme look like?

u0[{x_,y_}]:=Cos[x] Sin[y];
uInit=SetIC[{nx,ny,Count2,p},u0];
u=uInit;
Do[
u = LinearSolve[ANew, AOld.u ],
{12}];

ListContourPlot[Partition[u,nx+3]]
ListContourPlot[Partition[uInit,nx+3]]

ListPlot3D[Partition[u,nx+3]]
ListPlot3D[Partition[uInit,nx+3]]

?SetIC