Heat Equation Crank-Nicolson Scheme

MATH 330 Differential Equations II

Run the following to clear all variables:

Clear["Global`*"]

Parameters

Length of rod is len, thermal diffusivity coefficient is γ, time step is Δt, number of spatial intervals is n. 
Assign some values to these parameters here:

len=1;
γ=1.1/2; 
Δt = .01; 
n=50;

From the previous parameters we can compute Δx and μ:

Δx=(len/n);
μ=((γ Δt)/(((Δx))^2));

Initial conditions

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

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

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]=N[Table[ f[m*Δx],{m,1,n-1}]]

Boundary conditions

For now, we will use zero Dirichlet boundary conditions (that is, zero temperature at the endpoints). In finite difference notation, (Subscript[u, j,0])=0 and (Subscript[u, j,n])=0 for all time indexes j. We will see how these boundary conditions come into play below.

Matrices

The Crank-Nicolson scheme gives the following equation:

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

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^((j+1)))=B (u^((j))). We will compute (u^((j+1))) using the matrix equation (u^((j+1)))=(A^(-1))B (u^((j))).

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:

To build matrices A and B, we create (n - 1) by (n - 1) zero matrices and then use Do loops to set values along three diagonals. Replace the question marks below with the correct entries of each matrix.

A = ConstantArray[0, {n-1,n-1}];
Do[A[[j,j]]=??,{j,1,n-1}]; (* set entries on the main diagonal *)
Do[A[[j,j+1]]=??,{j,1,n-2}];  (* set entries just above the main diagonal *)
Do[A[[j+1,j]]=??,{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]]=??,{j,1,n-1}]; (* set entries on the main diagonal *)
Do[B[[j,j+1]]=??,{j,1,n-2}];  (* set entries just above the main diagonal *)
Do[B[[j+1,j]]=??,{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 j in u[j]. As mentioned before, we use the equation (u^((j+1)))=(A^(-1))B (u^((j))). The following loop does all the work:

Do[u[j+1]=Inverse[A].B.u[j],{j,0,100}]

Also append our vector u to include zero temperature at the boundaries x = 0 and x = 1:

Do[PrependTo[u[j],0],{j,0,100}]
Do[AppendTo[u[j],0],{j,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]

Plot the approximate solution

Since u[m] returns a list , we can use ListPlot to plot the solution at a particular time.  The option Joined draws a solid line instead of plotting discrete points, and the option DataRange scales the x-axis properly.

ListPlot[u[10],Joined->True,DataRange->{0,1},AxesLabel->{"x","temperature"}]

We can also make a bunch of plots and animate them:

plots = Table[ListPlot[u[j],Joined->True,DataRange->{0,1},AxesLabel->{"x","temperature"},PlotRange->{0,1}],{j,0,80}];
ListAnimate[plots]

We can get fancy by animating the plots.

ListAnimate[plots]

Questions

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

Experiment with different initial conditions. What do you observe?

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