Heat Equation: finite difference approximation via the Crank Nicolson Scheme

Run the following to clear all variables:

Clear["Global`*"]

Parameters

Length of rod is L, coefficient of the heat equation is k, time step is Δt, number of spatial intervals is n. We assign some values to these parameters here:

L=1;
k=1.1/2; 
Δt = .01; 
n=10;
Δx=L/n; 
s=k Δt / (Δx)^2;

Initial conditions

Let the initial temperature distribution be f(x)=sin(π log((x+1)e/2)).

f[x_]:=Sin[Pi Log[(x+1)E/2]];
Plot[f[x],{x,0,L}]

We will evaluate f(x) at intervals of Δx and store the values in U[0], with the zero denoting time zero. This will give the initial values for our approximate solution.

U[0]=Table[ f[j*Δx],{j,1,n-1}]

Boundary conditions

For now, we will use zero Dirichlet boundary conditions (that is, zero temperature at the endpoints). Note that in the second boundary condition, the input is n, not L. We don’t need to do anything with these boundary conditions yet.

Matrices

The Crank-Nicolson scheme gives the following equation:

-(s/2)(Subscript[u, j-1]^((m+1)))+(1+s)(Subscript[u, j]^((m+1)))-(s/2)(Subscript[u, j+1]^((m+1)))=(s/2)(Subscript[u, j-1]^((m)))+(1-s)(Subscript[u, j]^((m)))+(s/2)(Subscript[u, j+1]^((m)))

We build two matrices A and B that encode the coefficients. Matrix A encodes the coefficients on the left, and B encodes the coefficients on the right. Each of these matrices is tridiagonal. In matrix form, the equation above becomes A (u^((m+1)))=B (u^((m))). We will compute (u^((m+1))) using the matrix equation (u^((m+1)))=(A^(-1))B (u^((m))).

Since we have spatial points (Subscript[x, 0]) to (Subscript[x, n]), but the boundary points have fixed value zero, we have n-1 interior points to solve for at each time.

To build matrix A, we create an (n - 1) by (n - 1) zero matrix, and then use Do loops to set values along three diagonals:

A = ConstantArray[0, {n-1,n-1}];
Do[A[[j,j]]=(1+s),{j,1,n-1}]; (* set entries on the main diagonal *)
Do[A[[j,j+1]]=-s/2,{j,1,n-2}];  (* set entries just above the main diagonal *)
Do[A[[j+1,j]]=-s/2,{j,1,n-2}]; (* set entries just below the main diagonal *)
MatrixForm[A]

Similarly, build the matrix B:

B = ConstantArray[0, {n-1,n-1}];
Do[B[[j,j]]=(1-s),{j,1,n-1}]; (* set entries on the main diagonal *)
Do[B[[j,j+1]]=s/2,{j,1,n-2}];  (* set entries just above the main diagonal *)
Do[B[[j+1,j]]=s/2,{j,1,n-2}]; (* set entries just below the main diagonal *)
MatrixForm[B]

Compute the approximate solution

Next we can iterate and find the numerical solution.  Let's go from m=0 to m=100.  Since Δt=.01, this corresponds to t = 0 to t = 1. We store the solution values at time step m in U[m]. As mentioned before, we use the equation (u^((m+1)))=(A^(-1))B (u^((m))). The following loop does all the work:

Do[U[m+1]=Inverse[A].B.U[m],{m,0,100}]

Inspect some values of our approximate solution

Examine the solution at a particular time, such as m=10, which corresponds to t=0.1:

U[10]

We can append our vector U to include zero temperature at the boundaries x = 0 and x = 1:

Do[PrependTo[U[m],0],{m,0,100}]
Do[AppendTo[U[m],0],{m,0,100}]

Now check U[10] again, and compare with what we had above:

U[10]

Plot the approximate solution

Since U[m] returns a list , we can use ListPlot to plot the solution at a particular time.  The command Joined draws a solid line instead of plotting discrete points.

ListPlot[U[0]]

ListPlot[U[0],Joined->True]

A nice way of visualizing the solution is to make a table of plots.  Below we plot the curves at times t = 0, .2, .4, .6 .8.  These times correspond to m=0,20,40,60,80.

plots = Table[ListPlot[U[m],Joined->True,PlotRange->{0,1}],{m,0,80,20}]

Note that in the above, the plot is ‘shifted’ to start at 1 and end at 11.  We can adjust this with the DataRange modifier.

plots = Table[ListPlot[U[m],Joined->True,PlotRange->{0,1},DataRange->{0,1}],{m,0,80,20}]

We can get fancy by animating the plots.

ListAnimate[plots]

We can get even fancier by making more plots to smooth out our animation.

betterplots=Table[ListPlot[U[m],Joined->True,PlotRange->{0,1},DataRange->{0,1}],{m,0,80}];
ListAnimate[betterplots]

Questions

Try running the Crank-Nicholson scheme with different choices of Δx and Δt. Does the scheme always appear to be stable?

How would you incorporate nonhomogeneous time-dependent boundary conditions? For example: u(0,t)=a(t) and u(L,t)=b(t) for some functions a(t) and b(t)? Try it!