4620 Homework 6
Travis Beaulieu
Greg Mink

Counter Equation

TwoDCounterForDirchletNummanPeriodic[i_,j_]:=1 + i + (nx+3)*Mod[j,ny+1]

Interior and Periodic Points Equation

We set the bottom row of periodic points equal to the top row to eliminate the top row.
The equation for an interior point is:
	(Subscript[UNow, i,j])-(Subscript[UPrevious, i,j])=(1/2)((Δt)/((Δx)^2))((Subscript[UNow, i-1,j])-2(Subscript[UNow, i,j])+(Subscript[UNow, i+1,j]))+(1/2)((Δt)/((Δy)^2))((Subscript[UNow, i,j-1])-2(Subscript[UNow, i,j])+(Subscript[UNow, i,j+1]))+(1/2)((Δt)/((Δx)^2))((Subscript[UPrevious, i-1,j])-2(Subscript[UPrevious, i,j])+(Subscript[UPrevious, i+1,j]))+(1/2)((Δt)/((Δy)^2))((Subscript[UPrevious, i,j-1])-2(Subscript[UPrevious, i,j])+(Subscript[UPrevious, i,j+1])) + Δt (1/2)((Subscript[fPrevious, i,j])+(Subscript[fNeumannow, i,j]))

Where UNow is the 2D array of values, uNow is the vector you get by squishing UNow flat, and  samething for UPrevious and uPrevious.  Rearranging we get:
(1)	(Subscript[UNow, i,j])-(1/2)((Δt)/((Δx)^2))((Subscript[UNow, i-1,j])-2(Subscript[UNow, i,j])+(Subscript[UNow, i+1,j]))-(1/2)((Δt)/((Δy)^2))((Subscript[UNow, i,j-1])-2(Subscript[UNow, i,j])+(Subscript[UNow, i,j+1]))
(2)	(Subscript[UPrevious, i,j])+(1/2)((Δt)/((Δx)^2))((Subscript[UPrevious, i-1,j])-2(Subscript[UPrevious, i,j])+(Subscript[UPrevious, i+1,j]))+(1/2)((Δt)/((Δy)^2))((Subscript[UPrevious, i,j-1])-2(Subscript[UPrevious, i,j])+(Subscript[UPrevious, i,j+1])) + ((Δt)/2)((Subscript[fPrevious, i,j])+(Subscript[fNeumannow, i,j]))
Also:
{αx,αy}=(1/2){((Δt)/((Δx)^2)),((Δt)/((Δy)^2))}

InteriorPointMatrixBuild[{i_,j_,TwoDCounterForDirchletNummanPeriodic_},{αx_,αy_},{AN_,AO_}] := Module[
{k=TwoDCounterForDirchletNummanPeriodic[i,j],ANow=AN, APrevious=AO},
ANow[[k,k]]=1+2αx + 2αy;
ANow[[k,TwoDCounterForDirchletNummanPeriodic[i-1,j]]]=-αx;
ANow[[k,TwoDCounterForDirchletNummanPeriodic[i+1,j]]]=-αx;
ANow[[k,TwoDCounterForDirchletNummanPeriodic[i,j-1]]]=-αy;
ANow[[k,TwoDCounterForDirchletNummanPeriodic[i,j+1]]]=-αy;
APrevious[[k,k]]=1-2αx - 2αy;
APrevious[[k,TwoDCounterForDirchletNummanPeriodic[i-1,j]]]=αx;
APrevious[[k,TwoDCounterForDirchletNummanPeriodic[i+1,j]]]=αx;
APrevious[[k,TwoDCounterForDirchletNummanPeriodic[i,j-1]]]=αy;
APrevious[[k,TwoDCounterForDirchletNummanPeriodic[i,j+1]]]=αy;
{ANow,APrevious}
]

Dirichlet Points Equation

The equation for a Dirichlet point is: 
	(Subscript[UNow, i,j])= fDirchlet((Subscript[x, i]),(Subscript[y, j]),(Subscript[t, New])).

DirichletPointMatrixBuild[{i_,j_,TwoDCounterForDirchletNummanPeriodic_},{αx_,αy_},{AN_,AO_}] := Module[
{k=TwoDCounterForDirchletNummanPeriodic[i,j],ANow=AN, APrevious=AO},
ANow[[k,k]]=1;
{ANow,APrevious}
]

Ghost Points Equation

The equation for a Ghost point is:
	(Subscript[UNow, i,j])-(Subscript[UNow, i-2,j])=2Δx fNeumann((Subscript[x, i-1]),(Subscript[y, j]),(Subscript[t, New])).

GhostPointMatrixBuild[{i_,j_,TwoDCounterForDirchletNummanPeriodic_},{αx_,αy_},{AN_,AO_}] := Module[
{k=TwoDCounterForDirchletNummanPeriodic[i,j],ANow=AN, APrevious=AO},
ANow[[k,k]]=1.0;
ANow[[k,TwoDCounterForDirchletNummanPeriodic[i-2,j]]]=-1.0;
{ANow,APrevious}
]

Forcing Terms

Setting the iteration to the form:
	ANow.uNow=APrevious.uPrevious+(OverVector[f])
where (OverVector[f]) is defined in terms of our input fEquation, fDirchlet, and fNeumann depending on which sort of point we are at.

?Do

UpDatef[{fEquation_,fDirchlet_,fNeumann_},{Δx_,Δy_,Δt_,m_}, {nx_,ny_}]:=Module[
(* Something missing? *)
{i,j,k,n,p,f},
n=(nx+3)(ny+1);
p[i_,j_]={a+i Δx, c+j Δy}; (* physical {x,y} coordinates for i,j *)
count[i_,j_]= 1 + i + Mod[j,ny+1]*(nx+3); (* matrix location for i,j *)
(* intializing f *)
f=Table[0,n];
(* updating Dirichlet conditions *)
i=0;
Do[
k=count[i,j];
f[[k]]=fDirchlet[p[i,j],(m+1)*Δt],
{j,0,ny}
];
(* updating Neumann conditions *)
i=nx+2;
Do[
k=count[i,j];
f[[k]]=2*Δx*fNeumann[p[i-1,j],(m+1)*Δt],
{j,0,ny}
];
(* updating interior conditions *)
Do[
k=count[i,j];
f[[k]]=0.5*Δt(fEquation[p[i,j],m*Δt]+fEquation[p[i,j],(m+1)*Δt]),
{i,1,nx+1},
{j,0,ny}
];
f
]

Initial Conditions

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

Putting it all together

{{a,b},{c,d}}={{0,2},{3, 4}};
{nx, ny}={3, 5};
{Δx,Δy}={((b-a)/(nx+1)),((d-c)/(ny+1))};
p[i_,j_]:={a+i Δx,c + j Δy}
TwoDCounterForDirchletNummanPeriodic[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))};
ANow=SparseArray[{},{n,n}];
APrevious=SparseArray[{},{n,n}];

MatrixBuilder[{nx_,ny_},{αx_,αy_},{ANow_,APrevious_}]:=Module[
(* Something missing? *)
(* To create the interior points of A *)
{i,j,TwoDCounterForDirchletNummanPeriodic},
TwoDCounterForDirchletNummanPeriodic[i_,j_]= 1 + i + Mod[j,ny+1](nx+3);
Do[
{ANow,APrevious}=InteriorPointMatrixBuild[{i,j,TwoDCounterForDirchletNummanPeriodic},{αx,αy},{ANow,APrevious}],
{i,1,nx+1},
{j, 0, ny}
];
(* To create the dirichlet points of A *)
Do[
{ANow,APrevious}=DirichletPointMatrixBuild[{ 0,j,TwoDCounterForDirchletNummanPeriodic},{αx,αy},{ANow,APrevious}],
{j, 0, ny}
];
(* To create the ghost points of A *)
Do[
{ANow,APrevious}=GhostPointMatrixBuild[{nx+2,j,TwoDCounterForDirchletNummanPeriodic},{αx,αy},{ANow,APrevious}],
{j, 0, ny}
];
{ANow,APrevious}
]

{MatrixPlot[APrevious],MatrixPlot[ANow]}

{,}

u0[{x_,y_}]:=Cos[2. π y] (x^(3/2)+2);
uInitial=SetInitialCondition[{nx,ny,TwoDCounterForDirchletNummanPeriodic,p},u0];
u=uInitial;
Do[
u = LinearSolve[ANow, APrevious.u ],
{16}];

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

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

Uanalytical[x_,y_,t_]=Cos[2 π y] E^(-t) (x^3+2)
fEquation=Cos[2 π y] E^(-t) ((2 π)^2 (x^3+2)-(x^3+2)-(6 x+2))
fDirchlet=2 Cos[2 π y] E^(-t)
fNeumann=14 Cos[2 π y] E^(-t)
UpDatef[{fEquation_,fDirchlet_,fNeumann_}]