Heat Equation Explicit 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=10;

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

Mathematica stores exact values whenever possible. For faster computation we will use the N (“numeric”) function to store decimal values:

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.

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 step.

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

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

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-2μ) and the adjacent diagonals to μ. 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-2μ),{j,1,n-1}];
Do[A[[j,j+1]]=μ,{j,1, n-2}];
Do[A[[j+1,j]]=μ,{j,1,n-2}];
MatrixForm[A]

Compute the approximate solution

Next we can iterate and find the numerical solution.  Let's go from time index j=0 to time index j=100.  Since Δt=.01, this corresponds to time values from t=0 to t=1. We store the solution values at time step m in u[j]. 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,99}]

Inspect some values of our approximate solution

Examine the solution at a particular time index, such as j=10, which corresponds to time (Subscript[t, 10])=0.1:

u[10]

We can append the zero boundary values to our vector u like this:

Do[PrependTo[u[j],0],{j,0,100}]
Do[AppendTo[u[j],0],{j,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 option Joined draws a solid line instead of plotting discrete points. Try it:

ListPlot[u[0]]

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

Note that in the plot above, the list points are plotted at values 1,2,...,11 on the horizontal axis. We can use the DataRange option to fix this. We will also add axis labels.

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

A nice way of visualizing the solution is to make a table of plots.  Below we plot the curves at times t=0,0.2,0.4,0.6,0.8.  These times correspond to time indexes j=0,20,40,60,80.

plots = Table[ListPlot[u[j],Joined->True,DataRange->{0,1},AxesLabel->{"x","temperature"},PlotRange->{0,1}],{j,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.

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

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(t,0)=a(t) and u(t,len)=b(t) for some functions a(t) and b(t)? Try it!