CN Periodic BC 1 Spatial Dimension

Mp[n_,ν_]:= SparseArray[{
(* Periodic Solution Matrix for CN n is the number of interior pts*)
(* Includes the zeroth point but not the other end *)
Band[{1,1}]->1-2 ν,
Band[{1,2}]->ν,
Band[{2,1}]->ν,
{n+1,1}->ν,
{1,n+1}->ν
},
{n+1,n+1}]

CNSolver[uIn_List,{Δt_, nT_}] :=Module[
{n=Length[uIn]-1, Δx,ν,MPlus, MMinus,IntervalLength=2π,u},
u=uIn;
Δx=IntervalLength/(n+1);
ν=Δt/((Δx)^2);
MPlus  =Mp[n,   0.5 ν];
MMinus=Mp[n,-0.5ν];
(* Solve MMinus.(U^(n+1))= MPlus.(U^n) repeatedly storing solution in u *) 
Do[
u=LinearSolve[MMinus,MPlus.u],
{nT}];
u
]

uTest1[x_,t_]:=2(E^(-(3^2)t)) Sin[3x]+ 4 (E^(-(2^2)t))Cos[2x] ;
nx=96;
Δx =2π/(nx+1);
u0=Table[uTest1[x,0],{x,0.0,2π-Δx,Δx}] ;
{Δt, nt}={0.001,1000};
uEnd=CNSolver[u0,{Δt, nt}];
Show[
Plot[uTest1[x, nt*Δt],{x,0, 2π}],
ListPlot[uEnd,
DataRange->{0, 2π - Δx}]
]

a

Time Trials

n=(10^3);
SpA = Mp[n,0.1];
FA = Normal[SpA];
b = RandomReal[{-1,1},n+1];
Timing[xF=LinearSolve[FA,b];]
Timing[B=Inverse[FA];]
Do[B.b,{40}]//Timing
Do[xSp=LinearSolve[SpA,b],{40}]//Timing
LUDecomposition[SpA];//Timing
Norm[xSp-xF]

Clear[f]
f=LinearSolve[SpA];//Timing
z=Table[0,{i,500}];
Do[z[[i]]=f[i*b],{i,1,500}]//Timing

2 Spatial Dimensions - Periodic

First we need to work out how we are going to label our grid points.  We decided to use the simple wraparound recipe shown below.  e are not computing the two sides shown (we need to change this if we are doing non-peridic boundar conditions).

Clear[i,j,nx,ny,ax,ay,bx,by,x,y,count]
{nx,ny}={6,4};
{ax,bx}={0,π};
{ay,by}={0, π};
{Δx,Δy}= {bx-ax,by-ay}/{nx+1,ny+1};
x[i_]= ax+i Δx;
y[j_]=ay + j Δy;
count[i_,j_]= 1+i + (nx+1)*j;
Graphics[{
Table[
{Red,Tooltip[Point[{x[i],y[j]}],{i,j}]},
{i,0,nx+1},
{j,0,ny+1}
],
Table[
Text[count[i,j],{x[i],y[j]},{-1,-1}],
{i,0,nx},
{j,0,ny}
]
},
Frame->True,
FrameLabel->{"x","y"}
]

dim=(nx+1)(ny+1);
(* Seting up a blank matrix of the correct size *)
A=SparseArray[{},{dim,dim}];
(* Filling in the generic entries *)
Do[
eq=count[i,j];
A[[eq,eq]]=-4.0;
A[[eq,count[i+1,j]]]=1.0;
A[[eq,count[i-1,j]]]=1.0;
A[[eq,count[i,j+1]]]=1.0;
A[[eq,count[i,j-1]]]=1.0,
{i,1,nx-1},
{j,1,ny-1}
]
(* filling in the non-corner i=0 entries *)
Clear[i,j]
i=0;
Do[
eq=count[i,j];
A[[eq,eq]]=-4.0;
A[[eq,count[i+1,j]]]=1.0;
A[[eq,count[nx, j]]]=1.0; (* periodic boundary condition *)
A[[eq,count[i,j+1]]]=1.0;
A[[eq,count[i,j-1]]]=1.0,
{j,1,ny-1}
]
(* filling in the non-corner i=nx entries *)
Clear[i,j]
i=nx;
Do[
eq=count[i,j];
A[[eq,eq]]=-4.0;
A[[eq,count[0,j]]]=1.0; (* periodic boundary condition *)
A[[eq,count[i-1, j]]]=1.0;
A[[eq,count[i,j+1]]]=1.0;
A[[eq,count[i,j-1]]]=1.0,
{j,1,ny-1}
]
(* filling in the non-corner j=ny entries *)
Clear[i,j]
j=ny;
Do[
eq=count[i,j];
A[[eq,eq]]=-4.0;
A[[eq,count[i+1,j]]]=1.0;
A[[eq,count[i-1,j]]]=1.0;
A[[eq,count[i,0       ]]]=1.0; (* periodic boundary condition *)
A[[eq,count[i,j-1]]]=1.0,
{i,1,nx-1}]
(* filling in the non-corner j=0 entries *)
Clear[i,j]
j=0;
Do[
eq=count[i,j];
A[[eq,eq]]=-4.0;
A[[eq,count[i+1,j]]]=1.0;
A[[eq,count[i-1,j]]]=1.0;
A[[eq,count[i,j+1]]]=1.0;  
A[[eq,count[i,ny    ]]]=1.0, (* periodic boundary condition *)
{i,1,nx-1}]
(* filling in the four corners {0,0}, {nx,0}, {O, ny}, {nx, ny} *)
Clear[i,j]
{i,j}={0,0};
eq=count[i,j];
A[[eq,eq]]=-4.0;
A[[eq,count[i+1,j]]]=1.0;
A[[eq,count[nx,j]]]=1.0; (* periodic boundary condition *)
A[[eq,count[i,j+1]]]=1.0;  
A[[eq,count[i,ny    ]]]=1.0; (* periodic boundary condition *)
Clear[i,j]
{i,j}={nx,0};
eq=count[i,j];
A[[eq,eq]]=-4.0;
A[[eq,count[0,       j]]]=1.0; (* periodic boundary condition *)
A[[eq,count[i-1,j]]]=1.0; 
A[[eq,count[i,j+1]]]=1.0;  
A[[eq,count[i,ny    ]]]=1.0; (* periodic boundary condition *)
Clear[i,j]
{i,j}={nx,ny};
eq=count[i,j];
A[[eq,eq]]=-4.0;
A[[eq,count[0,       j]]]=1.0; (* periodic boundary condition *)
A[[eq,count[i-1,j]]]=1.0; 
A[[eq,count[i,0       ]]]=1.0; (* periodic boundary condition *)  
A[[eq,count[i,j-1]]]=1.0; 
Clear[i,j]
{i,j}={0,ny};
eq=count[i,j];
A[[eq,eq]]=-4.0;
A[[eq,count[i+1,       j]]]=1.0; 
A[[eq,count[nx,           j]]]=1.0;  (* periodic boundary condition *)
A[[eq,count[i,0       ]]]=1.0; (* periodic boundary condition *)  
A[[eq,count[i,j-1]]]=1.0; 


MatrixForm[A]
MatrixPlot[A]

Some sort of check!

{evals,evecs}=Eigensystem[Normal[A]];
evals[[-1]]
evecs[[-1]]

Looking at the eigenvectors

Manipulate[
ListPlot3D[Partition[evecs[[i]],{nx+1}],
PlotLabel->StringForm["i=`` and eval=``",i,evals[[i]]]
],
{i,Length[evecs],1,-1}]

Another way of doing the same thing!

If we think a bit harder about the edge pieces we can probably make it work automatically and save ourselves a lot of typing.  As Jack said think modulus.  The question is can we make it automatically put in the correct modified count.  The answer is "Of Course".  The key is to look back at the guide and realize that you want nx+1->0 and nx+1->0

?Mod

Table[
Mod[i,nx+1],
{i,0, 2 nx}]

dim=(nx+1)(ny+1);
(* Seting up a blank matrix of the correct size *)
AA=SparseArray[{},{dim,dim}];
(* Filling in the generic entries *)
Do[
eq=count[i,j];
AA[[eq,eq]]=-4.0;
AA[[eq,count[Mod[i+1,nx+1],Mod[j,       ny+1]]]]=1.0;
AA[[eq,count[Mod[i-1,nx+1],Mod[j,       ny+1]]]]=1.0;
AA[[eq,count[Mod[i,      nx+1], Mod[j+1,ny+1]]]]=1.0;
AA[[eq,count[Mod[i,      nx+1], Mod[j-1,ny+1]]]]=1.0,
{i,0,nx},
{j,0,ny}
]


MatrixForm[AA]
MatrixPlot[AA]

Norm[A-AA]

3 Spatial Dimensions - Periodic

Clear[i,j,k]
{nx,ny,nz}={16,17,18};
dim=(nx+1)(ny+1)(nz+1);
count[i_,j_,k_]:= 1+i + (nx+1)*j + (nx+1)(ny+1)*k;
(* Seting up a blank matrix of the correct size *)
AA=SparseArray[{},{dim,dim}];
(* Filling in the generic entries *)
Do[
eq=count[i,j,k];
AA[[eq,eq]]=-6.0;
AA[[eq,count[Mod[i+1,nx+1],Mod[j,       ny+1],Mod[k,         nz+1]]]]=1.0;
AA[[eq,count[Mod[i-1,nx+1],Mod[j,       ny+1],Mod[k,         nz+1]]]]=1.0;
AA[[eq,count[Mod[i,      nx+1], Mod[j+1,ny+1],Mod[k,         nz+1]]]]=1.0;
AA[[eq,count[Mod[i,      nx+1], Mod[j-1,ny+1],Mod[k,         nz+1]]]]=1.0;
AA[[eq,count[Mod[i,      nx+1], Mod[j       ,ny+1],Mod[k+1,  nz+1]]]]=1.0;
AA[[eq,count[Mod[i,      nx+1], Mod[j       ,ny+1],Mod[k-1,  nz+1]]]]=1.0,
{i,0,nx},
{j,0,ny},
{k,0, nz}
]

MatrixPlot[AA]

(nx+1)(ny+1)

Timing[
f=LinearSolve[AA]
]

b=RandomReal[{-1,1},Length[AA]];
b=b-Mean[b];
Timing[Do[f[i*b],{i,1,100}]]