Heat Equation: finite difference approximation via the Richardson 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;

Richardson Scheme

The Richardson scheme uses the following equation to step forward in time:

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

where s=((k Δt)/((Δx)^2)). Recall that (Subscript[u, j]^((m)))denotes the approximation at time step m (that is, time m Δt) and position index j (that is, position (Subscript[x, j])).

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}]

Since Richardson’s scheme requires values at two consecutive time steps to compute values at the next time step, the values in U[0] aren’t enough to get started. One workaround is simply to store these same values in U[-1]. (This might not be the best idea, but it should work. Can you think of a better solution?)

U[-1]=U[0]

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

Recall the equation for Richardson’s scheme:

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

We build a matrix A that encodes the coefficients in the equation. This will allow us to rewrite the equation in vector form as (u^((m+1)))=A (u^((m)))+(u^((m-1))), which will be useful for computation.

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]]=-4s,{j,1,n-1}]; (* set entries on the main diagonal *)
Do[A[[j,j+1]]=2s,{j,1,n-2}];  (* set entries just above the main diagonal *)
Do[A[[j+1,j]]=2s,{j,1,n-2}]; (* set entries just below the main diagonal *)
MatrixForm[A]

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 (u^((m)))+(u^((m-1))). The following loop does all the work:

Do[U[m+1]=A.U[m]+U[m-1],{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

What happens if you run Richardson’s scheme with different choices of Δx or Δt? Do you observe any behavior that appears stable?

In order to compute the first iteration of Richardson’s scheme, the code above stored the values of U[0] in U[-1]. Is this a reasonable thing to do? Can you think of a better solution?