4 An explicit numerical method for the heat equation

The approximations used above for approximating partial derivatives can now be applied in order to derive a numerical method for solving the heat conduction problem

$\begin{array}{ccccc}\hfill u_t& \hfill =\hfill & \alpha u_{xx}\hfill & \hfill \hfill & \left(0<x<\ell ,t>0\right)\hfill \\ \hfill u\left(0,t\right)& \hfill =\hfill & 0\hfill & \hfill \hfill & \left(t>0\right)\hfill \\ \hfill u\left(\ell ,t\right)& \hfill =\hfill & 0\hfill & \hfill \hfill & \left(t>0\right)\hfill \\ \hfill u\left(x,0\right)& \hfill =\hfill & f\left(x\right)\hfill & \hfill \hfill & \left(0<x<\ell \right).\hfill \end{array}$

In order to specify the numerical method we choose values for $\delta t$ and $\delta x$ and use these in approximations of the two derivatives in the partial differential equation. It is convenient to divide the interval $0<x<\ell$ into equally spaced subintervals so, in effect, we choose a whole number $J$ so that $\delta x=\frac{\ell }{J}$ .

Figure 3

The diagram above shows the independent variables $x$ and $t$ at which we seek the function $u$ . The numerical solution we shall find is a sequence of numbers which approximate $u$ at a sequence of $\left(x,t\right)$ points.

The notation we use is that

$\begin{array}{ccc}\hfill u_j^n\hfill & \hfill \approx \hfill & \hfill \underset{⏟}{u\left(j\delta x,n\delta t\right)}\hfill \\ \hfill \uparrow \hfill & \hfill \hfill & \hfill \uparrow \hfill \\ \hfill \text{numerical}\hfill & \hfill \hfill & \hfill \text{exact (i.e., unknown) solution}\hfill \\ \hfill \text{approximation}\hfill & \hfill \hfill & \hfill \text{evaluated at }x=j\times \delta x\text{, }t=n\times \delta t\hfill \end{array}$

The idea is that the subscript $j$ counts how many “steps" to the right we have taken from the origin and the superscript $n$ counts how many time-steps (up, on the diagram) we have taken. To say this another way

For example, consider the point on Figure 3 which is highlighted with a small square. This point is two steps to the right of the origin (so that $j=2$ ) and five steps up (so that $n=5$ ). The exact solution evaluated at this point is $u\left(2\delta x,5\delta t\right)$ and our numerical approximation to that value is $u_2^5$ .

Combining this new notation with the familiar idea for approximating derivatives we obtain the following approximation to the PDE

$\frac{u_j^{n+1}-u_j^n}{\delta t}=\alpha \frac{u_{j-1}^n-2u_j^n+u_{j+1}^n}{{\left(\delta x\right)}^2}$

To simplify (the appearance of) the numerical method we define a new quantity    $r=\frac{\alpha \delta t}{{\left(\delta x\right)}^2}$   so that our numerical procedure can be written

$u_j^{n+1}=u_j^n+r\left(u_{j-1}^n-2u_j^n+u_{j+1}^n\right)=r{u}_{j-1}^n+\left(1-2r\right)u_j^n+r{u}_{j+1}^n$

This equation defines a numerical “stencil" which allows us to find one of the values at the $n+1$ time level in terms of values at the previous level, $n$ . In Figure 4 we envisage terms on the right-hand side of the above equation leading towards a result equal to the left-hand side, and the arrows therefore point towards the point at which $u_j^{n+1}$ approximates $u$ .

Figure 4

At the stage of the process depicted above, the solid circles represent points in the $\left(x,t\right)$ plane where we have already found our numerical approximation. The unfilled circle is the point for which the new approximation $u_j^{n+1}$ is being found.

4.1 Implementation

The initial condition gives $u$ at $t=0$ , and this information can be used to find

$u_0^0,u_1^0,u_2^0,\dots ,u_J^0$

that is, the numerical solution at all the selected $x$ values and at $t=0$ . In general

$u_j^0=f\left(j\times \delta x\right)=f_j$

where $f_j$ is a shorthand notation for $f\left(j\times \delta x\right)$ .

Then we use the boundary conditions and numerical method

$u_j^{n+1}=u_j^n+r\left(u_{j-1}^n-2u_j^n+u_{j+1}^n\right)$

(with $n=0$ ) to work out $u_j^1$ for $j=0,1,2,\dots ,J$ . (This completes the first time-step.)

The time-stepping procedure is then used repeatedly to find $u_j^{n+1}$ in terms of the $u_j^n$ , which are known either from the last time-step or (at the beginning) from the initial condition.

The time-stepping procedure is summarised in the following Key Point.

The following is a concrete example of the time-stepping procedure.

Example 15

The temperature $u\left(x,t\right)$ of a metal bar of length $\ell =2$ at a distance $x$ from one end and at time $t$ is modelled by the partial differential equation

$u_t=\alpha u_{xx}\left(0<x<\ell ,t>0\right)$

It is given that the metal has diffusivity $\alpha =4$ , that the two ends of the bar are kept at temperature $u=0$ and that the initial temperature distribution is

$u\left(x,0\right)=f\left(x\right)=x\left(\ell -x\right)$

Use the explicit difference scheme with $\delta x=0.5$ and $\delta t=0.01$ to approximate $u\left(x,t\right)$ at $t=\delta t$ and $t=2\delta t$ .

Solution

In this case $r=\alpha \delta t∕{\left(\delta x\right)}^2=0.16$ so that the numerical method can be written

$u_j^{n+1}=u_j^n+0.16\left(u_{j-1}^n-2u_j^n+u_{j+1}^n\right)=0.68u_j^n+0.16\left(u_{j-1}^n+u_{j+1}^n\right)$

We now find $u_j^0$

$\begin{array}{cccc}\hfill u_0^0& \hfill =\hfill & 0\hfill & \text{from the left-hand boundary condition}\hfill \\ \hfill u_1^0& \hfill =\hfill & f\left(\delta x\right)=0.75\hfill & \text{from the initial condition}\hfill \\ \hfill u_2^0& \hfill =\hfill & f\left(2\delta x\right)=1\hfill & \text{from the initial condition}\hfill \\ \hfill u_3^0& \hfill =\hfill & f\left(3\delta x\right)=0.75\hfill & \text{from the initial condition}\hfill \\ \hfill u_4^0& \hfill =\hfill & 0\hfill & \text{from the boundary condition at the right hand end}\hfill \end{array}$

The first time-step will find $u_j^1$ , but first we note that $u_0^1=u_4^1=0$ from the two boundary conditions. Now

$\begin{array}{ccc}\hfill u_1^1& \hfill =\hfill & 0.68u_1^0+0.16\left(u_0^0+u_2^0\right)=0.68\times 0.75+0.16\left(0+1\right)=0.670\hfill \\ \hfill u_2^1& \hfill =\hfill & 0.68u_2^0+0.16\left(u_1^0+u_3^0\right)=0.68\times 1+0.16\left(0.75+0.75\right)=0.920\hfill \\ \hfill u_3^1& \hfill =\hfill & 0.68u_3^0+0.16\left(u_2^0+u_4^0\right)=0.68\times 0.75+0.16\left(1+0\right)=0.670\hfill \\ \hfill \end{array}$

The second time-step will find $u_j^2$ , but first we note that $u_0^2=u_4^2=0$ from the two boundary conditions. Now

$\begin{array}{ccc}\hfill u_1^2& \hfill =\hfill & 0.68u_1^1+0.16\left(u_0^1+u_2^1\right)=0.68\times 0.67+0.16\left(0+0.92\right)=0.603\hfill \\ \hfill u_2^2& \hfill =\hfill & 0.68u_2^1+0.16\left(u_1^1+u_3^1\right)=0.68\times 0.92+0.16\left(0.67+0.67\right)=0.84\hfill \\ \hfill u_3^2& \hfill =\hfill & 0.68u_3^1+0.16\left(u_2^1+u_4^1\right)=0.68\times 0.67+0.16\left(0.92+0\right)=0.603\hfill \\ \hfill \end{array}$

(Quantities have been rounded to three decimal places here.)

Figure 5 plots the numerical solutions found in the example above. The initial condition is shown as circles. Results of the first time-step appear as squares and the second time-step is shown as stars. The line joining the values we found are not part of the numerical solution and are included only as an aid to clarity.

Figure 5

Notice how the numerical results are behaving as they should. The temperature decreases slightly at each time-step.

Task!

The temperature $u\left(x,t\right)$ of a metal bar of length $\ell =2$ at a distance $x$ from one end and at time $t$ is modelled by the partial differential equation

$u_t=\alpha u_{xx}\left(0<x<\ell ,t>0\right)$

It is given that the metal has diffusivity $\alpha =2.25$ , that the two ends of the bar are kept at temperature $u=0$ and that the initial temperature distribution is

$u\left(x,0\right)=f\left(x\right)=sin\left(\pi x∕\ell \right)$

Use the explicit difference scheme with $\delta x=0.5$ and $\delta t=0.05$ to approximate $u\left(x,t\right)$ at $t=\delta t$ and $t=2\delta t$ .

Answer

Answer