$$ \newcommand{\D}{\displaystyle} \renewcommand{\eqref}[1]{Eq.~(\ref{#1})} $$




7.5.4 Generalized von Neumann stability analysis

We have referred to the stability condition \( |G|\leq1 \) as the von Neumann's strong stability condition. The reason for this designation is that if \( |G|\leq 1 \) is fulfilled, then the amplitude of the solution can not increase. However, there are many physical problems for which amplitude grows in time (in a bounded manner). A simple example is that of heat conduction with a source such as: $$ \begin{equation} \tag{7.73} \frac{\partial T}{\partial t}= \alpha\frac{\partial^2 T}{\partial x^2}+bT,\ b= const.,\ t < t_{max} \end{equation} $$

\( b < 0 \) represents a heat sink, while \( b > 0 \) is said to be a heat source. In the first case we can apply the von Neumann strong criterion, whereas in the later case it is necessary to allow for \( |G|>1 \).

A particular solution of (7.73) is given by: $$ \begin{equation} \tag{7.74} T(x,t)=e^{bt}\cdot e^{-\alpha\beta^2\cdot t}\cos(\beta x)=e^{(b-\alpha\beta^2)\cdot t}\cos(\beta x) \end{equation} $$

Now, let us use (7.74) to determine an analytical growth factor, see (7.39): $$ \begin{equation} \tag{7.75} G_a=\frac{T(x_j,\ t_{n+1})}{T(x_j,\ t_n)}=\exp\left[(b-\alpha\beta^2)\cdot t_{n+1}-(b-\alpha\beta^2)\cdot t_n\right]=e^{b\Delta t}\cdot e^{-\alpha\beta^2\Delta t} \end{equation} $$

We see that the term \( e^{b\cdot \Delta t} \) causes increases of amplitude for positive \( b \).

A series expansion for small \( \Delta t \): $$ \begin{equation} \tag{7.76} e^{b\Delta t}=1+b\cdot \Delta t+ \frac{b^2}{2}(\Delta t)^2+\dots \end{equation} $$

If we use the FTCS-scheme we obtain: $$ \begin{equation} \tag{7.77} T^{n+1}_{j}=D(T^{n}_{j+1}+T^{n}_{j-1})+(1-2D)T^{n}_{j}+\Delta t bT^{n}_{j} \end{equation} $$

with $$ \begin{equation} \tag{7.78} D=\alpha \frac{\Delta t}{(\Delta x)^2} \end{equation} $$

If we use the PC-criterion we find that the sum of coefficients is equal to \( 1+b\cdot \Delta t \), which implies that this criterion can only be used for \( b < 0 \).

On the other hand, using the von Neumann's method yields (7.77): $$ \begin{equation} \tag{7.79} G=1-4D\sin^2\Big(\frac{\delta}{2}\Big)+\Delta t\cdot b \end{equation} $$

By setting \( D=\frac{1}{2} \), we have: $$ \begin{equation*} |G|\leq \left|1-2\sin^2\Big(\frac{\delta}{2}\Big)\right|+|\Delta t\cdot b|\leq 1+\Delta t\cdot b,\ b > 0 \end{equation*} $$

If we compare this expression with (7.76), we see that the agreement between the analytical and the numerical growth factor for this case.

We can now introduce the generalized von Neumann's stability condition: $$ \begin{equation} \tag{7.80} |G|\leq 1+K\cdot \Delta t \end{equation} $$ with \( K \) being a positive constant.

This means that we allow the amplitude to grow exponentially for \( t < t_{max} \). In this case we can use the strong stability condition if we reason this way: since the source term in (7.73) doe not contain any derivative term, we can ignore this term in the stability analysis. We have the same kind of problem in section 1.6 in the Numeriske Beregninger, where stiff ordinary differential equations are discussed (See for example equation (1.6.8) in section 1.6.1 of the Numeriske Beregninger).

For a growing amplitude , we must decrease the step length of the independent variable if we are to achieve a prescribed accuracy. For a decreasing amplitude, however, we must keep it under a maximum step length to get a stable calculation.

Let us look at another example using the FTCS-scheme. $$ \begin{equation} \tag{7.81} \frac{\partial u}{\partial t}=\alpha \frac{\partial^2 u}{\partial x^2}+a_0\frac{\partial u}{\partial x},\ \alpha > 0 \end{equation} $$

This equation is called the advection-diffusion equation and is a parabolic equation according to the classification scheme discussed in section 4.3 of the Numeriske Beregninger.

Using the FTCS scheme for (7.81) with central differences for \( \dfrac{\partial u}{\partial x} \): $$ \begin{equation*} u_j^{n+1}=u^n_{j}+D\cdot(u^n_{j+1}-2u^n_{j}+u^n_{j-1})+a_0\frac{\Delta t}{2\Delta x}(u^n_{j+1}-u^n_{j-1}),\ D=\alpha \frac{\Delta t}{(\Delta x)^2} \end{equation*} $$

Using von Neumann's method we get that: $$ \begin{equation*} G=1-4D\sin^2\Big(\frac{\delta}{2}\Big) +i\cdot a_0 \frac{\Delta t}{\Delta x}\sin(\delta) \end{equation*} $$

which further provides: $$ \begin{equation*} |G|^2=\left(1-4D\sin^2\Big(\frac{\delta}{2}\Big)\right)^2+\left(a_0\frac{\Delta t}{\Delta x}\sin(\delta)\right)^2=\left(1-4D\sin^2\Big(\frac{\delta}{2}\Big)\right)^2+\frac{a_0^2\cdot D}{\alpha}\cdot \Delta t\sin^2(\delta ) \end{equation*} $$

Setting now \( D=\frac{1}{2} \): $$ \begin{equation*} |G|=\sqrt{\left(1-2\sin^2\Big(\frac{\delta}{2}\Big)\right)^2+\frac{a_0^2}{2\alpha}\cdot\Delta t\cdot \sin^2(\delta)} \leq 1+\frac{a_0^2}{2\alpha}\cdot \Delta t \end{equation*} $$

Here we have used the inequality \( \sqrt{x^2+y^2}\leq |x|+|y| \)

With \( K=\frac{a^2_0}{2\alpha} \), we see that the generalized condition (7.80) is met.

The advection-diffusion equation is discussed in detail in section6.10 of the Numeriske Beregninger.