5 Stability of the simple explicit scheme

The purpose of the time-stepping scheme is to approximate $u\left(x,t\right)$ at later and later times $t$ . It is clear that the larger we take the time step $\delta t$ , the fewer steps will be necessary to reach a particular time $t$ . One constraint on the size of $\delta t$ is that we know from our earlier look at difference methods that derivative approximations are most accurate when small increments are used. However, as we will see in the next couple of pages, a far more telling constraint on the size of $\delta t$ arises on consideration of stability . We begin with an Example.

Example 16

The temperature $u\left(x,t\right)$ of a metal bar of length $\ell =1$ 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 =1$ , 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.25$ and $\delta t=0.075$ 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=1.2$ so that the numerical scheme can be written

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

The first stage is to use the given data to find $u_j^0$

$\begin{array}{ccccc}\hfill u_0^0& \hfill =\hfill & 0\hfill & \hfill & \text{from the boundary condition}\hfill \\ \hfill u_1^0& \hfill =\hfill & f\left(\delta x\right)=f\left(0.25\right)=0.188\hfill & \hfill & \text{from the initial condition}\hfill \\ \hfill u_2^0& \hfill =\hfill & f\left(2\delta x\right)=f\left(0.5\right)=0.25\hfill & \hfill & \text{from the initial condition}\hfill \\ \hfill u_3^0& \hfill =\hfill & f\left(3\delta x\right)=f\left(0.75\right)=0.188\hfill & \hfill & \text{from the initial condition}\hfill \\ \hfill u_4^0& \hfill =\hfill & 0\hfill & \hfill & \text{from the boundary condition}\hfill \end{array}$

The first time-step will find $u_j^1$ . First we note that the boundary condition implies that $u_0^1=u_4^1=0$ .

$\begin{array}{ccc}\hfill u_1^1& \hfill =\hfill & -1.4u_1^0+1.2\left(u_0^0+u_2^0\right)=-1.4\times 0.19+1.2\left(0+0.25\right)=0.038\hfill \\ \hfill u_2^1& \hfill =\hfill & -1.4u_2^0+1.2\left(u_1^0+u_3^0\right)=-1.4\times 0.25+1.2\left(0.188+0.188\right)=0.1\hfill \\ \hfill u_3^1& \hfill =\hfill & -1.4u_3^0+1.2\left(u_2^0+u_4^0\right)=-1.4\times 0.19+1.2\left(0.25+0\right)=0.038\hfill \\ \hfill \end{array}$

The second time-step will find $u_j^2$ . First we note that the boundary condition implies that $u_0^2=u_4^2=0$ . Now

$\begin{array}{ccc}\hfill u_1^2& \hfill =\hfill & -1.4u_1^1+1.2\left(u_0^1+u_2^1\right)=-1.4\times 0.04+1.2\left(0+0.1\right)=0.067\hfill \\ \hfill u_2^2& \hfill =\hfill & -1.4u_2^1+1.2\left(u_1^1+u_3^1\right)=-1.4\times 0.1+1.2\left(0.038+0.038\right)=-0.05\hfill \\ \hfill u_3^2& \hfill =\hfill & -1.4u_3^1+1.2\left(u_2^1+u_4^1\right)=-1.4\times 0.04+1.2\left(0.1+0\right)=0.067\hfill \\ \hfill \end{array}$

Figure 6 shows the results found in Example 16.

Figure 6

Something has gone wrong here. And it only gets worse in subsequent time-steps. After 9 time-steps the numerical solution approximating $u\left(x,t\right)$ at $t=9\delta t$ is

(to 4 decimal places). This is an example of instability . A part of the numerical solution wants to keep growing and growing in a way that is not a part of the engineering application being modelled.

There are many different definitions of (in)stability, and they often depend on the specific application in mind. For the heat conduction problem under discussion here, the following definition is sufficient.

(Of course, there are applications where the principal quantity of interest does grow with time, and in these cases other definitions of stability are appropriate.)

The main stability result for the explicit scheme is proved in many textbooks on the subject, but for this Workbook it is sufficient to simply state it.

5.1 Why is the stability constraint a problem?

In the above account it has been stated that the stability constraint is a severe restriction on the time-step $\delta t$ . Here we discuss why this is the case.

For sake of argument let us take an example where $\alpha =1$ and choose $\delta x=\frac{1}{10}$ . The stability requirement insists that we must choose

$\delta t\le \frac{1}{2}\delta x^2=\frac{1}{200},$

which is much smaller than $\delta x$ . If we require an even smoother approximation in the $x$ direction we could halve $\delta x$ taking it to be equal to $\frac{1}{20}$ . It is now necessary that

$\delta t\le \frac{1}{2}\delta x^2=\frac{1}{800}.$

Decreasing $\delta x$ by a factor of 2 causes $\delta t$ to decrease by a factor of 4. The problem is that the upper bound on $\delta t$ involves the square of $\delta x$ , which is likely to be very small.

The following method overcomes the requirement of tiny time-steps.