Heat Equation: finite difference approximation via the explicit scheme

The explicit scheme is also called the Forward-Time Central Space (FTCS) 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}]

Mathematica stores exact values whenever possible. To see decimal approximations, use the N (“numeric”) function:

U[0]=N[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 will see how these boundary conditions come into play below.

Matrix

Next we build the matrix A that will allow us to easily compute the approximation values at each successive time step. 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.

We first create an (n - 1) by (n - 1) zero matrix:

A = ConstantArray[0, {n-1,n-1}]

We will use a forward difference approximation for the time variable, and centered difference for the spatial variable. As we have seen, we need to set the main diagonal to (1-2s) and the adjacent diagonals to s. We call the entries of A using the syntax A[[i,j]], and we use three Do loops to set the values:

Do[A[[j,j]]=(1-2s),{j,1,n-1}];
Do[A[[j,j+1]]=s,{j,1, n-2}];
Do[A[[j+1,j]]=s,{j,1,n-2}];
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]. The following loop does all the work, but it doesn’t print any output! (Note that the dot denotes matrix/vector multiplication.)

Do[U[m+1]=A.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 print 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. Try it:

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

Does our solution look like what you expected? Is something wrong? If so, how would you fix it? Experiment!

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!