Making Sparse Matrices

n=3;
m=3;
A=SparseArray[
Join[
Table[ 
{i,i}->-4,
{i, m*n}],
Table[{i,i+1}->1,{i,m*n-1}],
Table[{i+1,i}->1,{i,m*n-1}],
Table[{i+m,i}->1,{i,m*(n-1)}],
Table[{i,i+m}->1,{i,m*(n-1)}]
]
]
MatrixForm[A]
Det[A]
MatrixPlot[A]

Ist Discretization of the 1D Heat Equation

Consider (Subscript[u, t])= κ (Subscript[u, xx])+f  on a<x<b and 0<t<T with Δx=((b-a)/mx)and Δt=(T/my) subject to the initial condition u[x,0]=g[x] and the Dirichlet boundary conditions u[a,t]=u[b,t]=0.  Writing U[[i,n]]≈u[(Subscript[x, i]),(Subscript[t, n])] where (Subscript[x, i])=b+i Δx and (Subscript[t, n])=0+n Δt and using a forward difference in t and a centered difference in x at u[(Subscript[x, i]),(Subscript[t, n])] gives the finite difference approximation

u[(Subscript[x, i]),(Subscript[t, n+1])]-u[(Subscript[x, i]),(Subscript[t, n+1])] = ((κ Δt)/((Δx)^2))(u[(Subscript[x, i+1]),(Subscript[t, n])]-2u[(Subscript[x, i]),(Subscript[t, n])]+u[(Subscript[x, i-1]),(Subscript[t, n])])+Δt f[(Subscript[x, i]),(Subscript[t, n])]

with initial conditions

u[(Subscript[x, i]), 0]=g[(Subscript[x, i])]

and boundary conditions

u[(Subscript[x, 0]), (Subscript[t, j])]=u[(Subscript[x, mx]), (Subscript[t, n])]=0.

In terms of our matrix U (and α=((κ Δt)/((Δx)^2))) this is

U[[i,n+1]]-U[[i,n]] = α(U[[i+1, n]]-U[[i, n]]+U[[i+1, n]]) + Δt F[[i,n]]

where F[[i,n]]=f[(Subscript[x, i]),(Subscript[t, n])] with boundary conditions

U[[0,n]]=U[[mx,n]]=0

and initial conditions

U[[i,0]]=g[(Subscript[x, i])].

Note: (4) can be solved explicitly for the unknown value of U[[i,n+1]] and we just march up through the j values.

Note, we need to decide if we start counting from zero or one in the matrix.  As an example, lets compute the approximate solution to (1) with κ=1, f[x,t]=(x^2)t and initial condition g[x]=0 on the domain with a=0, b=1,and T=10 with mx=3 and mt=5.

First the set up.

f[x_,t_]=(x^2)t;
κ=1;
{a,b}={0.0,1.0};
{mx,mt}={10,2001};
T=10.0;
{Δx,Δt}={((b-a)/mx),(T/mt)};
α = κ ((Δt)/((Δx)^2));
F=Table[
f[x, t], 
{x,a,b,Δx},
{t,0,T,Δt}];
MatrixPlot[F]

Now the run.

U=Table[
0, 
{x,a,b,Δx},
{t,0,T,Δt}];
Do[
Do[
U[[i,j+1]] = U[[i,j]]+α ( U[[i+1,j]] - 2 U[[i,j]] + U[[i-1,j]]) +Δt F[[i,j]],
{i,2,mx}],
{j,1,mt}
]
MatrixPlot[U]
ListPlot3D[U]

Here is a plot of the solution at time t=T i.e. All the x values and the last time value.

ListPlot[U[[All,-1]]]

2nd Discretization of the 1D Heat Equation

Consider (Subscript[u, t])= κ (Subscript[u, xx])+f  on a<x<b and 0<t<T with Δx=((b-a)/mx)and Δt=(T/mt) subject to the initial condition u[x,0]=g[x] and the Dirichlet boundary conditions u[a,t]=u[b,t]=0.  Writing U[[i,n]]≈u[(Subscript[x, i]),(Subscript[t, n])] where (Subscript[x, i])=b+i Δx and (Subscript[t, n])=0+n Δt and using a backward difference in t and a centered difference in x at u[(Subscript[x, i]),(Subscript[t, n+1])] gives the finite difference approximation

u[(Subscript[x, i]),(Subscript[t, n+1])]-u[(Subscript[x, i]),(Subscript[t, n])] = ((κ Δt)/((Δx)^2))(u[(Subscript[x, i+1]),(Subscript[t, n+1])]-2u[(Subscript[x, i]),(Subscript[t, n+1])]+u[(Subscript[x, i-1]),(Subscript[t, n+1])])+Δt f[(Subscript[x, i]),(Subscript[t, n+1])]

with initial conditions

u[(Subscript[x, i]), 0]=g[(Subscript[x, i])]

and boundary conditions

u[(Subscript[x, 0]), (Subscript[t, j])]=u[(Subscript[x, mx]), (Subscript[t, j])]=0.

In terms of our matrix U (and α=((κ Δt)/((Δx)^2))) this is

U[[i,n+1]]-U[[i,n]] = α(U[[i+1, n+1]]-U[[i, n+1]]+U[[i+1, n+1]]) + Δt F[[i,n+1]]

where F[[i,n]]=f[(Subscript[x, i]),(Subscript[t, n])] with boundary conditions

U[[0,j]]=U[[n,0]]=0

and initial conditions

U[[i,0]]=g[(Subscript[x, i])].

Note: (10) can be written as a linear equation for the unknown vector of values U[[All,j+1]] in terms of the known values U[[All,j]].  It looks like this

A.U[[All,j+1]]=Δt F[[All,j+1]] +U[[All,j]]

where the matrix A=Id - α M with M our favorite 1,-2,1,… stencil matrix from the centered difference operation.  We need to think again about how to enforce the boundary conditions here!

f[x_,t_]=(x^2)t;
κ=1;
{a,b}={0.0,1.0};
{n,m}={23,124};
T=10.0;
{Δx,Δt}={((b-a)/n),(T/m)};
α = κ ((Δt)/((Δx)^2));
F=Table[
f[x, t], 
{x,a,b,Δx},
{t,0,T,Δt}];
MatrixPlot[F]

Now the run. First we set up the matrix for the second difference operator.

A=Table[0,
{n+1},
{n+1}
];
Do[A[[i,i]]=-2,{i,2,n}];
Do[A[[i,i-1]]=1;A[[i,i+1]]=1,{i,2,n}];
MatrixForm[A]

Now we set up the upDate matrix including the Identity and the boundary conditions by zeroing the boundary heat supplies

Id = IdentityMatrix[n+1];
UpDateMatrix=Id-α A;
F[[All,1]]=0*F[[All,1]];
F[[All,-1]]=0*F[[All,-1]];
MatrixForm[UpDateMatrix]
MatrixForm[F]

Now the run!

U=Table[
0.0, 
{x,a,b,Δx},
{t,0,T,Δt}
];

Do[
U[[All,j+1]] =LinearSolve[ UpDateMatrix,U[[All,j]]+Δt F[[All,j]]],
{j,1,m-1}
]
MatrixPlot[U]
ListPlot3D[U]

Accuracy: Local Truncation Error

Truncation Error For Explicit Scheme

The book defines the truncation error as simply the taylor series of what you get when you substitute the solution u[x,t] into the numerical scheme. Here is the result for the explicit scheme.

LTEExplicit=CoefficientList[
Normal[
Series[
((u[x,t+Δt]-u[x,t])/(Δt))-α ((u[x-Δx,t]-2u[x,t]+u[x+Δx,t])/((Δx)^2)),
{Δt,0, 3},
{Δx,0,3}
]
],
{Δt,Δx}
];
Grid[
LTEExplicit,
Frame->All
]

Note that at least one of these terms simplifies because of the heat equation.

Grid[
LTEExplicit/.{
D[u[x,t],t]-> α D[u[x,t],x,x]
},
Frame->All
]

Truncation Error For Implicit Scheme

The book defines the truncation error as simply the taylor series of what you get when you substitute the solution u[x,t] into the numerical scheme. Here is the result for the explicit scheme.

LTEImplicit=CoefficientList[
Normal[
Series[
((u[x,t]-u[x,t-Δt])/(Δt))-α ((u[x-Δx,t]-2u[x,t]+u[x+Δx,t])/((Δx)^2)),
{Δt,0, 3},
{Δx,0,3}
]
],
{Δt,Δx}
];
Grid[
LTEImplicit,
Frame->All
]

Note that at least one of these terms simplifies because of the heat equation.

Grid[
LTEImplicit/.{
D[u[x,t],t]-> α D[u[x,t],x,x]
},
Frame->All
]

Truncation Error For θ Scheme

The θ scheme simply uses a weighted average (with weight θ ) of the implicit and explicit schemes:  θ=1 is fully implicit and θ=0 is fully explicit.

LTEθ=CoefficientList[
Normal[
Series[
((u[x,t+Δt]-u[x,t])/(Δt))-(1-θ)α ((u[x-Δx,t]-2u[x,t]+u[x+Δx,t])/((Δx)^2))
-
θ α ((u[x-Δx,t+Δt]-2u[x,t+Δt]+u[x+Δx,t+Δt])/((Δx)^2))
,
{Δt,0, 3},
{Δx,0,3}
]
],
{Δt,Δx}
];
Grid[
Simplify[LTEθ],
Frame->All
]

Note that at least one of these terms simplifies because of the heat equation.

Grid[
LTEθ/.{
D[u[x,t],t]-> α D[u[x,t],{x,2}]
},
Frame->All
]

Grid[
Simplify[
LTEθ/.{
D[u[x,t],t]-> α D[u[x,t],{x,2}],
D[u[x,t],{t,2}]-> α D[u[x,t],{x,2},{t,1}]
}
],
Frame->All
]

Comment on the choice θ=1/2.

Stability: Error Propagation

Our analysis of stability is motivated by the existence of special solutions of the heat equation with no heat source.  It is easy to check that

(Subscript[u, n])(x,t)=(e^(- α (n^2) t))sin(n x)

is a solution of the heat equation (Subscript[u, t])=(Subscript[u, xx]) with Dirichlet boundary conditions on the interval 0≤x≤π.  Note, the interval of length π is just for convenience.  Summing a bunch of these Fourier solutions gives us a fairly general solution of the form.

u(x,t)=(∑_(n=1)^∞)(Subscript[b, n])(e^(- α (n^2) t))sin(n x)

If we want to do the periodic problem (again of length π) then we need to add in cosines and a constant (which we could think of as cos(0 x) if we wanted to) to get the form

u(x,t)=(Subscript[a, 0])+(∑_(n=1)^∞)(Subscript[a, n])(e^(- α (n^2) t))sin(n x)+(Subscript[b, n])(e^(- α (n^2) t))sin(n x).

which (since (e^(i θ))=cos(θ)+i sin(θ)) we can combine in the complex form

u(x,t)=(∑_(n=0)^∞)(Subscript[a, n])(e^(- α (n^2) t))(e^(i n x)).

The coefficient (Subscript[a, n]) and (Subscript[b, n]) are Fourier coefficents of a periodic function. It turns out that "any" periodic function (for instance a periodic initial condition f) can be represented as a Fourier Series

f(x)=(∑_(n=0)^∞)(Subscript[a, n])(e^(i n x))

and these values of (Subscript[a, k]) immediately define solution of our heat equation.  The big question is what does this have to do with our difference schemes?  Well our difference schemes have analogous solutions

(Subscript[U, j]^n)=(λ^n)(e^(i k (j Δx)))

Substituting this form into our explicit scheme (not to bad if we note that

(Subscript[U, j]^(n+1))=λ (Subscript[U, j]^n)

(Subscript[U, j+1]^n)- 2(Subscript[U, j]^n) +(Subscript[U, j-1]^n) =(Subscript[U, j]^n)((e^(i k Δx))-2+(e^(-i k Δx)))

gives a solution with

λ = 1-4μ ((sin(1/2 k Δx))^2)