Organization

We need organization.

1D Problem

Our one spatial dimension {a≤x≤b} problem is easy to organize:
  Physical:		i=0,1,…,nx   with (Subscript[x, i])=a+i*Δx where Δx=(b-a)/nx
  Linear Algebra:	k=1,…,nx+1 with
  Connection:		k=i+1   
  			or  i=k-1
  Interior Pts:		i=1,…,nx-1
  			or k=2,…,nx
  Boundary Pts:	i=0 (left) and i=nx (right)
  			or k=1(left)  and k=nx+1(right)
  Ghost Pts:		i=-1(left)  and i=nx+1(right)
  (if needed)		or k=0 (left) and k=nx+2(right)
  
  Interior Eq:		1≤i≤nx-1 
  			or 2≤k≤nx
  (explicit)		(Subscript[u, i]^(n+1))=(Subscript[u, i]^n)+νx((Subscript[u, i-1]^n)-2 (Subscript[u, i]^n)+(Subscript[u, i+1]^n))+Δt (Subscript[f, i]^n)
  (implicit)		(Subscript[u, i]^(n+1))=(Subscript[u, i]^n)+νx((Subscript[u, i-1]^(n+1))-2 (Subscript[u, i]^(n+1))+(Subscript[u, i+1]^(n+1)))+Δt (Subscript[f, i]^(n+1))
  (θ)			(Subscript[u, i]^(n+1))=(Subscript[u, i]^n)+ νx (θ((Subscript[u, i-1]^n)-2 (Subscript[u, i]^n)+(Subscript[u, i+1]^n))+(1-θ)((Subscript[u, i-1]^(n+1))-2 (Subscript[u, i]^(n+1))+(Subscript[u, i+1]^(n+1))))+Δt (Subscript[f, i]^(n+θ))
  (Crank-Nicholson)	(Subscript[u, i]^(n+1))=(Subscript[u, i]^n)+ νx ((1/2)((Subscript[u, i-1]^n)-2 (Subscript[u, i]^n)+(Subscript[u, i+1]^n))+(1/2)θ((Subscript[u, i-1]^(n+1))-2 (Subscript[u, i]^(n+1))+(Subscript[u, i+1]^(n+1))))Δt (Subscript[f, i]^(n+1/2))
  
   Boundary Eq:	i=0,nx
   			or k=1,nx+1	
  
  (Dirichlet)		(Subscript[u, 0]^(n+1))=(Subscript[Data, 0]^(n+1))  and/or (Subscript[u, nx]^(n+1))=(Subscript[Data, nx]^(n+1))
  u=Data
  
  (Neumann)		Interior Eq at i=0 with the ghost point eliminated using
  (Subscript[u, x])=NData		 (((Subscript[u, 0]^m)-(Subscript[u, -1]^m))/(Δx))=(Subscript[NData, 0]^m) for m=n,n+1, or n+1/2
  			(((Subscript[u, nx+1]^m)-(Subscript[u, nx]^m))/(Δx))=(Subscript[NData, nx]^m) for m=n,n+1, or n+1/2
  
  (Robin)		Interior Eq at i=0 with the ghost point eliminated using
  ux+αu=RData	 (((Subscript[u, 0]^m)-(Subscript[u, -1]^m))/(Δx))+α (Subscript[u, 0]^m)=(Subscript[RData, 0]^m) for m=n,n+1, or n+1/2
  			 (((Subscript[u, nx+1]^m)-(Subscript[u, nx]^m))/(Δx))+α (Subscript[u, nx]^m)=(Subscript[RData, nx]^m) for m=n,n+1, or n+1/2
  			 
  (Periodic)		Interior Eq at i=0 with the ghost points eliminated using warp around.

To build the matrices we can use the linear algebra index k or the physical index i. The only difference is the start and stop points. When specifying a vector or matrix entry I need to use the linear algebra counting i.  There are several ways to make the matrices for the 1D problem.

Messy (something to plan from)

nx = 4;
AInt=Table[0,{nx+1},{nx+1}];
Do[
k=i+1;
AInt[[k,k-1]]=1.0; AInt[[k,k]]=-2.0; AInt[[k,k+1]]=1.0,
{i,1,nx-1}
]
MatrixForm[AInt]

nx = 21;
{a,b}={0.0,3.14159};
u0[x_]= Sin[2x];
f[x_,t_]= Sin[3x] Sin[t];
κ=0.1;
Δx=((b-a)/nx);
Δt=0.1;
νx= κ((Δt)/((Δx)^2))
AInt=SparseArray[{},{nx+1,nx+1}];
F=Table[0,{nx+1}];
U0=Table[0,{nx+1}];
Do[
k=i+1; (* physical to LA index*)
AInt[[k,k-1]]=1.0; AInt[[k,k]]=-2.0; AInt[[k,k+1]]=1.0;U0[[k]]=u0[a + i Δx], (* Setting IC  *)
{i,1,nx-1}
]
Id=IdentityMatrix[nx+1]; Id[[1,1]]=0; Id[[nx+1,nx+1]]=0;

NIter=33;
Show[
ListLinePlot[U0,
DataRange->{a,b},
PlotRange->{-1,3}],
U=U0;
Do[
F = Table[f[a+i*Δx,(n-1) Δt],{i,0, nx}];
U =(Id + νx AInt).U + Δt*F,
{n,1,NIter}];
ListLinePlot[U,
DataRange->{a,b},
PlotStyle->Red],
PlotLabel->Δt*NIter
]

NIter=35;
Show[
ListLinePlot[U0,
DataRange->{a,b},
PlotRange->{-1,3}],
U=U0;
Do[
F = Table[f[a+i*Δx,(n-1) Δt],{i,0, nx}];
U =(Id + νx AInt).U + Δt*F,
{n,1,NIter}];
ListLinePlot[U,
DataRange->{a,b},
PlotStyle->Red],
PlotLabel->Δt*NIter
]

nx = 21;
{a,b}={0.0,3.14159};
u0[x_]= Sin[2x]+1 + 2 (x/π);
f[x_,t_]= 0 Cos[3x] Sin[t];
κ=0.1;
Δx=((b-a)/nx);
Δt=0.1;
νx= κ((Δt)/((Δx)^2))
AInt=SparseArray[{},{nx+1,nx+1}];
F=Table[0,{nx+1}];
U0=Table[0,{nx+1}];
Do[
k=i+1; (* physical to LA index*)
AInt[[k,k-1]]=1.0; AInt[[k,k]]=-2.0; AInt[[k,k+1]]=1.0;U0[[k]]=u0[a + i Δx], (* Setting IC  *)
{i,1,nx-1}
]
U0[[1]]=u0[a];
U0[[nx+1]]=u0[b];
Id=IdentityMatrix[nx+1]; Id[[1,1]]=0; Id[[nx+1,nx+1]]=0;

NIter=350;
Show[
ListLinePlot[U0,
DataRange->{a,b},
PlotRange->{-1,3}],
U=U0;
Do[
F =Δt* Table[f[a+i*Δx,(n-1) Δt],{i,0, nx}];
F[[1]]=1.0;F[[nx+1]]=3.0;
U =(Id + νx AInt).U + F,
{n,1,NIter}];
ListLinePlot[U,
DataRange->{a,b},
PlotStyle->Red],
PlotLabel->Δt*NIter
]

F[[1]]
F[[1]]

Slightly Organized

So I think I will make a couple of functions with obvious names that make rows of the matrix.  I need to make a function for the equations at internal points, the equations at Robin (Neumann is a special case of Robin) boundary points, and periodic points.

For right now I will assume I have the following defined:
	{a, b}, nx, {Δx,Δt}, νx, f, g,α and
	i is the physical spatial counter and n is the physical temporal counter, while
	k is the Linear Algebra counter.
We want to build the sparse array AExp∈(ℝ^((nx+1)×(nx+1))) and vector F∈(ℝ^(nx+1)) to implement the explicit iteration
	(u^new)=(Id + νx A) (u^old)+F
to approximate the solution of 
	(Subscript[u, t])=κ (Subscript[u, xx])+f
of a≤x≤b and t>0 with the intial condition
	u[x,0]=u0[x];
and either periodic conditions 
(P)	u[a-,t]=u[b+,t]
or one of the following boundary conditions at x=a
(D)	u[a, t] = g[a,t] 
(R/N)	(Subscript[u, x])[a, t]+α[a,t] u[a,t] = g[a,t] 
and one of the  boundary conditions at x=ab
(D)	u[b, t] = g[a,t] 
(R/N)	(Subscript[u, x])[b, t]+αb[b,t]u[b,t] = g[b,t].
Note:  α=0 reduces the Robin conditions to a Neumann condition.

We assume that the matrix A and vectors b and U0 are initialied to zero and add the appropriate entries.

InternalEq:= (
k=i+1;
A[[k,k-1]]=  1.0;
A[[k,k      ]]=-2.0;
A[[k,k+1]]=   1.0;
)
InternalRHS:=(
k=i+1;
F[[k]]=f[a+i*Δx,n*Δt];
)
DirichletBoundaryEq:=(
k=i+1;
Switch[ i,
0,   mult=  1; x=a,
nx, mult=-1; x=b
,_, Print["Error"]
];
A[[k,k]]=-1/νx;
)
DirichletRHS:=(
k=i+1;
Switch[ i,
0,   mult=  1; x=a,
nx, mult=-1; x=b
,_, Print["Error"]
];
F[[k]]=g[x,t];
)
RobinBoundaryEq:=(
k=i+1;
Switch[ i,
0,   mult=  1; x=a,
nx, mult=-1; x=b
,_, Print["Error"]
];
A[[k,k]]=-2.0+(1 + mult α[x] Δx);
A[[k,k+mult]]= 1.0;
)
RobinRHS:=(
k=i+1;
Switch[ i,
0,   mult=   1; x=a,
nx, mult=-1; x =b,
_, Print["Error"]
];
F[[k]]=f[x,n*Δt]- mult  Δx g[x,n*Δt];
)
PeriodicBoundaryEq:=(
k=i+1;
Switch[ i,
0, x=a,
nx,x=b,
_, Print["Error"]
];
A[[k,k]]=-2.0;
A[[k,2]]=  1.0;
A[[k,nx]]=1.0;
)
PeriodicRHS:=(
k=i+1;
Switch[ i,
0,   x=a,
nx, x =b,
_, Print["Error"]
];
F[[k]]=f[x,n*Δt];
)
IC:=(
k=i+1;
U0[[k]]=u0[a+i*Δx];
)

The stuff up above makes it much cleaner to set up and run a problem

Periodic Boundary Conditions

nx=12;
{a,b}={0.0,3.14159};
{Δx,Δt}={((b-a)/nx),0.01};
κ=0.5;
f[x_,t_]= Sin[x] Cos[t];
u0[x_]= Sin[x];
νx=κ((Δt)/((Δx)^2))
(* End of Domain intial conditon and forcing function specification *)
(* Intializations *)
A= SparseArray[{},{nx+1,nx+1}];
Id = IdentityMatrix[nx+1];
F = Table[0,{nx+1}];
U=F; U0=F;
(* Building matrices and RHS *)
n=0;
Do[
InternalEq; InternalRHS; IC,
{i,1,nx-1,1}];
Do[
PeriodicBoundaryEq; PeriodicRHS; IC,
{i,0,nx,nx}]
(* Start Run *)
U=U0;
nMax=4;
Do[
U=(Id + νx A).U + F; 
Do[InternalRHS,{i,1,nx-1}]; Do[PeriodicRHS,{i,0,nx,nx}],
{n,1,nMax}];
ListPlot[{U0,U},
DataRange->{a,b}
]

nMax=4;
Do[
U=(Id + νx A).U + F; 
Do[InternalRHS,{i,1,nx-1}]; Do[PeriodicRHS,{i,0,nx,nx}],
{n,1,nMax}];
ListPlot[{U0,U},
DataRange->{a,b}
]

nMax=54;
Do[
U=(Id + νx A).U + F; 
Do[InternalRHS,{i,1,nx-1}]; Do[PeriodicRHS,{i,0,nx,nx}],
{n,1,nMax}];
ListPlot[{U0,U},
DataRange->{a,b}]

Robin Boundary Condition at x=a and Dirichlet BC at x =b

nx=12;
{a,b}={0.0,3.14159};
{Δx,Δt}={((b-a)/nx),0.01};
κ=0.5;
f[x_,t_]= Sin[x] Cos[t];
g[a,t_]= 1.4;
g[b,t_]= 2.3;
α[a]=1.2;
u0[x_]= Sin[x];
νx=κ((Δt)/((Δx)^2))
(* End of Domain intial conditon and forcing function specification *)
(* Intializations *)
A= SparseArray[{},{nx+1,nx+1}];
Id = IdentityMatrix[nx+1];
F = Table[0,{nx+1}];
U=F; U0=F;
(* Building matrices and RHS *)
n=0;
Do[
InternalEq; InternalRHS; IC,
{i,1,nx-1,1}];
i=0;
RobinBoundaryEq; RobinRHS;
i=nx;
DirichletBoundaryEq; DirichletRHS;
(* Start Run *)
U=U0;

nMax=4;
Do[
U=(Id + νx A).U + F; 
Do[InternalRHS,{i,1,nx-1}]; 
i=0;
RobinBoundaryEq; RobinRHS;
i=nx;
DirichletBoundaryEq; DirichletRHS,
{n,1,nMax}];
ListPlot[{U0,U},
DataRange->{a,b}
]

U

More Organization

The huge pile of glabal variables is dangerous!  We should make functions so that we can isolate variables effectively.

Almost Organized

So we are almost done.  There are several problems with what we have done:  
1) We can easily mess ourselves up by changing a variable we do not intend to. 
2) We can not see the woods for the trees.
3) We have repetitive code blocks
I want to write four functions:
1) The first will build the matrix we need.   
2) The  second will build the RHS vector(s) we need. 
3) The third will build the IC we need.
4) The fourth will run the solver for a specified number of iterations.

First making matrices.

Basically all the boundary conditions do is modify an underlying matrix. 
The Periodic, Dirichlet and Robin matries are modifications of the underlying matrix.  The Dirichlet boundary condition needs to know the νx,  the Robin boundary condition needs to know the νx and the local coefficient α.

Now it is a pain to make a function for Dirichlet/Dirichlet, Robin/Dirichlet, Dirichlet/Robin, and Robin/Robin. We are going to combine them into one general boundary condition
	α (Subscript[u, x])+β u=g ⟶  α (((Subscript[u, 0])-(Subscript[u, -1]))/(Δx))+β (Subscript[u, 0])=g  and/or α (((Subscript[u, nx+2])-(Subscript[u, nx+1]))/(Δx))+β (Subscript[u, nx+1])=g
which gives the ghost points as 
	(Subscript[u, -1])=-((g Δx)/α) +((α+β Δx )/α)(Subscript[u, 1+nx])   and (Subscript[u, nx+2])=((g Δx)/α) +((α-β Δx )/α)(Subscript[u, 1+nx])
which we can substitute back into the internal equation to eleiminate the ghost point.  α=0 should give a Dirichlet condition if you look at this the formulas we have no problem letting α->0 but actually plugging in α=0 will break things.  It is traditional to simply use a very small value of α to enforce approximate Dirichlet conditions.  If you think about it this is what we were doing when we put a large value on the diagonals in HW1.

computation

Clear[α,β,i,Δx,nx,k]
fred0=Solve[α (((Subscript[u, 0])-(Subscript[u, -1]))/(Δx))+β (Subscript[u, 0])==g, (Subscript[u, -1])][[1]]//Simplify
Collect[k((((Subscript[u, -1])-2 (Subscript[u, 0])+(Subscript[u, 1])))/((Δx)^2))/.fred0, {(Subscript[u, 0]),(Subscript[u, 1])}]
fredn = Solve[α (((Subscript[u, nx+2])-(Subscript[u, nx+1]))/(Δx))+β (Subscript[u, nx+1])==g, (Subscript[u, nx+2])][[1]]//Simplify
Collect[k ((((Subscript[u, nx])-2 (Subscript[u, nx+1])+(Subscript[u, nx+2])))/((Δx)^2))/.fredn ,{(Subscript[u, nx]),(Subscript[u, nx+1])}]

Code

Clear[InternalA, PeriodicA,RNDA]
InternalA[nx_]:= Module[
{i,A=SparseArray[{},{nx+1,nx+1}]},
Do[
k=i+1;
A[[k,k-1]]=  1.0;A[[k,k      ]]=-2.0;A[[k,k+1]]=   1.0,
{i,1,nx-1}];
A
]
PeriodicA[nx_]:=Module[
{A,m},
A=InternalA[nx];
A[[nx+1,nx+1]]=-2.0; 
A[[nx+1,2]]=1.0;
A[[nx+1,nx]]=1.0;
A[[1,1]]=-2.0; 
A[[1,2]]=1.0;
A[[1,nx]]=1.0;
A
]
RNDA[nx_, {Δx_,{α0_,β0_,g0_},{αn_,βn_,gn_}}]:=Module[
{A=InternalA[nx],k},
k=1;
A[[k,k+1]]=1.0;A[[k,k]]=-2.0 +((α0+β0*Δx)/(α0)) ;
k=nx+1;
A[[k,k-1]]=1.0;A[[k,k]]=-2.0 + ((αn-βn*Δx)/(αn)) ;
A
]

These are not as efficient as they could be but it is much more important that we can test them until we are sure that they are correct.

Now the RHS

The RHS is also basically a modification of a simple evaluation

rhs[f_,t_,{a_,b_,Δx_}]:=Module[
{x},
Table[f[x,t],{x,a,b,Δx}]
]
Periodicrhs[f_,t_,{a_,b_,Δx_}]:= rhs[f,t,{a,b,Δx}]
RNDrhs[f_,t_,{a_,b_,Δx_},{κ_,{α0_,β0_,g0_},{αn_,βn_,gn_}}]:=Module[
{F=rhs[f,t,{a,b,Δx}]},
F[[   1]]-=((κ g0 )/(α0 Δx));
F[[-1]]+= ((κ gn )/(αn Δx));
F
]

Test: Periodic

Here is the periodic problem using the new functions

nx=36;
{a,b}={0.0,3.14159};
{Δx,Δt}={((b-a)/nx),0.001};
κ=0.5;
uSol[x_,t_]=Cos[2x]+3+Sin[2 x] Cos[x t];
f[x_,t_]=D[uSol[x,t],t]-κ D[uSol[x,t],{x,2}];
u0[x_]= uSol[x,0];
νx=κ((Δt)/((Δx)^2))
(* Building matrices and RHS and IC *)
A= PeriodicA[nx];
Id = IdentityMatrix[nx+1];
U0 =Table[u0[x],{x,a,b,Δx}];
TMax=1.2;
nMax=Floor[TMax/Δt]
(* Start Run *)
U=U0;
nMax=400;
Do[
F=Periodicrhs[f,n*Δt,{a,b,Δx}];
U=(Id + νx A).U +Δt F,
{n,0, nMax}] 
Show[
Plot[
uSol[x, nMax*Δt],
{x, a,b}],
ListPlot[{U},
DataRange->{a,b}
],
PlotRange->All
]

Test: Dirichlet/Neumann Homogneous

Here is the periodic problem using the new functions

nx=120;
{a,b}={0.0,3.14159};
{Δx,Δt}={((b-a)/nx),0.01};
κ=0.5;
f[x_,t_]= Sin[x] Cos[t];
u0[x_]= Sin[x];
νx=κ((Δt)/((Δx)^2))
(* Building matrices and RHS and IC *)
{{α0,β0,g0},{αn,βn,gn}}={{1,0,0},{1.0*(10^(-12)),1,0}}
A= RNDA[nx,{Δx,{α0,β0,g0},{αn,βn,gn}}];
Id = IdentityMatrix[nx+1];
U0 =Table[u0[x],{x,a,b,Δx}]
n=0;
(* Start Run *)
U=U0;
nMax=4;
Do[
F=RNDrhs[f,n*Δt,{a,b,Δx},{κ,{α0,β0,g0},{αn,βn,gn}}];
U=(Id + νx A).U + F,
{n,0, nMax}] 
ListPlot[{U0,U},
DataRange->{a,b}
]

Real Test

u[x_,t_]= 3 x+ ((Sin[3 x])^2)Sin[x- t]
nx=54;
{a,b}={0.0,3.141592653589793};
{Δx,Δt}={((b-a)/nx),0.001};
κ=0.5;
f[x_,t_]= D[u[x,t],t]-κ D[u[x,t],{x,2}];
u0[x_]= u[x,0];
Chop[D[u[x,t],x]/.{{x->a},{x->b}}]
νx=κ((Δt)/((Δx)^2))
(* Building matrices and RHS and IC *)
{{α0,β0,g0},{αn,βn,gn}}={{1,0,3},{1,0,3}}
A= RNDA[nx,{Δx,{α0,β0,g0},{αn,βn,gn}}];
Id = IdentityMatrix[nx+1];
U0 =Table[u0[x],{x,a,b,Δx}];
n=0;
(* Start Run *)
U=U0;
nMax=1000;
Do[
F=RNDrhs[f,n*Δt,{a,b,Δx},{κ,{α0,β0,g0},{αn,βn,gn}}];
U=(Id + νx A).U + Δt F,
{n,0, nMax}] 
Show[
Plot[ u[x,nMax*Δt],{x,a,b}],
ListPlot[{U,U0},
DataRange->{a,b},
PlotStyle->{Red,Green}
],
PlotRange->All
]

nMax*Δt