10.4: Stability

Another consideration for numerical schemes for the heat equation is the stability of the scheme. In implementing the finite difference scheme, \[u_{m, j+1}=u_{m, j}+\alpha\left[u_{m+1, j}-2 u_{m, j}+u_{m-1, j}\right]\label{eq:1}\] with \(\alpha=k \Delta t /(\Delta x)^{2}\), one finds that the solution goes crazy when \(\alpha\) is too big. In other words, if you try to push the individual time steps too far into the future, then something goes haywire. We can determine the onset of instability by looking at the solution of this equation for \(u_{m, j}\). [Note: We changed index \(i\) to \(m\) to avoid confusion later in this section.]

The scheme is actually what is called a partial difference equation for \(u_{m, j}\). We could write it in terms of difference, such as \(u_{m+1, j}-u_{m, j}\) and \(u_{m, j+1}-u_{m, j}\). The furthest apart the time steps are are one unit and the spatial points are two units apart. We can see this in the stencils in Figure 10.6. So, this is a second order partial difference equation similar to the idea that the heat equation is a second order partial differential equation. The heat equation can be solved using the method of separation of variables. The difference scheme can also be solved in a similar fashion. We will show how this can lead to product solutions.

We begin by assuming that \(u_{m j}=X_{m} T_{j}\), a product of functions of the indices \(m\) and \(j\). Inserting this guess into the finite difference equation, we have \[\begin{align} u_{m, j+1} &=u_{m, j}+\alpha\left[u_{m+1, j}-2 u_{m, j}+u_{m-1, j}\right]\nonumber \\ X_{m} T_{j+1} &=X_{m} T_{j}+\alpha\left[X_{m+1}-2 X_{m}+X_{m-1}\right] T_{j},\nonumber \\ \frac{T_{j+1}}{T_{j}} &=\frac{\alpha X_{m+1}+(1-2 \alpha) X_{m}+\alpha X_{m-1}}{X_{m}} .\label{eq:2} \end{align}\]

Noting that we have a function of \(j\) equal to a function of \(m\), then we can set each of these to a constant, \(\lambda\). Then, we obtain two ordinary difference equations: \[\begin{align} T_{j+1} &=\lambda T_{j},\label{eq:3} \\ \alpha X_{m+1}+(1-2 \alpha) X_{m}+\alpha X_{m-1} &=\lambda X_{m} .\label{eq:4} \end{align}\]

The first equation is a simple first order difference equation and can be solved by iteration: \[\begin{align} T_{j+1} &=\lambda T_{j},\nonumber \\ &=\lambda\left(\lambda T_{j-1}\right)=\lambda^{2} T_{j-1},\nonumber \\ &=\lambda^{3} T_{j-2},\nonumber \\ &=\lambda^{j+1} T_{0},\label{eq:5} \end{align}\]

The second difference equation can be solved by making a guess in the same spirit as solving a second order constant coefficient differential equation.Namely, let \(X_{m}=\xi^{m}\) for some number \(\xi\). This gives \[\begin{align} \alpha X_{m+1}+(1-2 \alpha) X_{m}+\alpha X_{m-1} &=\lambda X_{m},\nonumber \\ \xi^{m-1}\left[\alpha \xi^{2}+(1-2 \alpha) \xi+\alpha\right] &=\lambda \xi^{m}\nonumber \\ \alpha \xi^{2}+(1-2 \alpha-\lambda) \xi+\alpha &=0 .\label{eq:6} \end{align}\]

This is an equation foe \(\xi\) in terms of \(\alpha\) and \(\lambda\). Due to the boundary conditions, we expect to have oscillatory solutions. So, we can guess that \(\xi=|\xi| e^{i \theta}\), where \(i\) here is the imaginary unit. We assume that \(|\xi|=1\), and thus \(\xi=e^{i \theta}\) and \(X_{m}=\xi^{m}=e^{i m \theta}\). Since \(x_{m}=m \Delta x\), we have \(X_{m}=e^{i x_{m} \theta / \Delta x}\). We define \(\beta=\) theta / \(\Delta\), to give \(X_{m}=e^{i \beta x_{m}}\) and \(\xi=e^{i \beta \Delta x}\).

Inserting this value for \(\xi\) into the quadratic equation for \(\xi\), we have \[\begin{align} 0 &=\alpha \xi^{2}+(1-2 \alpha-\lambda) \xi+\alpha\nonumber \\ &=\alpha e^{2 i \beta \Delta x}+(1-2 \alpha-\lambda) e^{i \beta \Delta x}+\alpha\nonumber \\ &=e^{i \beta \Delta x}\left[\alpha\left(e^{i \beta \Delta x}+e^{-i \beta \Delta x}\right)+(1-2 \alpha-\lambda)\right]\nonumber \\ &=e^{i \beta \Delta x}[2 \alpha \cos (\beta \Delta x)+(1-2 \alpha-\lambda)]\nonumber \\ \lambda &=2 \alpha \cos (\beta \Delta x)+1-2 \alpha .\label{eq:7} \end{align}\]

So, we have found that \[u_{m j}=X_{m} T_{j}=\lambda^{m}\left(a \cos \alpha x_{m}+b \sin \alpha x_{m}\right), \quad a^{2}+b^{2}=h_{0}^{2},\nonumber \] and \[\lambda=2 \alpha \cos (\beta \Delta x)+1-2 \alpha\nonumber \] For the solution to remain bounded, or stable, we need \(|\lambda| \leq 1\).

Therefore, we have the inequality \[-1 \leq 2 \alpha \cos (\beta \Delta x)+1-2 \alpha \leq 1 .\nonumber \] Since \(\cos (\beta \Delta x) \leq 1\), the upper bound is obviously satisfied. Since \(-1 \leq\) \(\cos (\beta \Delta x)\), the lower bound is satisfied for \(-1 \leq-2 s+1-2 s\), or \(s \leq \frac{1}{2}\). Therefore, the stability criterion is satisfied when \[\alpha=k \frac{\Delta t}{\Delta x^{2}} \leq \frac{1}{2} .\nonumber \]