Heuristic Analysis
Table of Contents
Finite-difference equations may have rapidly growing and oscillating solutions that in no way resemble the solutions expected from the partial differential equations they are meant to approximate. Such solutions are said to exhibit computational instability. Clearly, it is desirable to avoid these numerical disasters. For linear difference equations with constant coefficients, computational stability can be determined using a Fourier method pioneered by von Neumann (see the article in this series “Computational Stability.” Unfortunately, most equations of physical interest are either nonlinear, or have non-constant coefficients, or both.
Heuristic Analysis Methods
In this article a simple heuristic analysis method is described for investigating the computational stability of such finite-difference equations. An important by-product of this type of analysis is that it often suggests simple ways to eliminate the instabilities and at the same time increase the accuracy of the approximations.
The approach described here is called “heuristic” because it is not rigorous or complete, but it often works and can provide a great deal of useful information. Reference [1] is the original publication describing the heuristic stability method from which much of this article has been taken.
Heuristic analysis is based on the rather simple idea of reducing a finite-difference equation back to a partial differential equation by expanding each of its terms in a Taylor series and keeping only terms to a certain order in the expansion. This expansion is in powers of the space and time increments, which are assumed to be small to begin with.
Certainly such an expansion must, to lowest order, reproduce the original partial differential equation, otherwise, it would not be a good approximation. Oftentimes this requirement is referred to as the “consistency” of the approximation. Terms beyond the lowest order in the expansion are referred to as truncation errors.
The basic concept of a heuristic analysis is that the Taylor-expanded equation is a more accurate representation of the difference equation than the original partial differential equation. Even keeping only a few truncation error terms should result in a partial differential equation that is more closely related to the difference equation. With this in mind, the following discussion will show that an examination of the truncated equation can sometimes reveal properties shared with the difference equation such as stability problems, necessary initial conditions and/or serious inaccuracies.
To begin, we consider the same linear partial differential equation that was discussed in the first article on stability: Computational Stability.
Linear Equation Example
The equation for one-dimensional advection-diffusion of a variable u(x,t) is
(1) $latex \displaystyle \frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=\nu \frac{{{\partial }^{2}}u}{\partial {{x}^{2}}}.$
The convection velocity c and the diffusion coefficient ν are assumed to be constants. Solutions of this equation are known to be bounded and otherwise well-behaved.
What will be shown here is that the stability of a simple finite-difference approximation to Eq. 1 can be determined from an examination of the truncations errors resulting from a Taylor series expansion of a the difference equation. Not only does this process reveal that there are two basic types of instability, but we shall be able to make a direct comparison between the heuristic method and the von Neumann type of Fourier analysis carried out in Computational Stability. This comparison provides a useful rule-of-thumb for which truncation error terms to keep and which to eliminate from the Taylor expansion in order to evaluate the difference equation’s stability.
The simple, explicit finite-difference equation approximating Eq. 1 discussed in Computational Stability is
(2) $latex \displaystyle \frac{u_{j}^{n+1}-u_{j}^{n}}{\delta t}=-\frac{c}{2\delta x}\left( u_{j+1}^{n}-u_{j-1}^{n} \right)+\frac{\nu }{\delta {{x}^{2}}}\left( u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n} \right)$
where, e.g., ujn denotes u(jδx,nδt). This is called a forward-in-time approximation that allows all j location values to be computed at time step n+1, provided all the j values at step n are known. In other words, the difference equation requires one initial condition to start things off, just as the original partial differential equation also requires a single initial condition because it only involves a single time derivative.
It may be observed that difference equation, Eq. 2, has the property that each space and time location (jδx,nδt) will affect points at time step n+1 at locations j-1, j and j+1. That is, point (jδx,nδt) has a region of influence at later time bounded by lines having slopes ±δx/δt in x-t space. These are similar to characteristic lines along which signals can propagate. For example, the original equation, Eq. 1, has a characteristic line with slope c along which a disturbance advects. In the discrete equation, however, the characteristic lines are not physical characteristics but computational ones defining the region where the difference equation changes data values resulting from a change in value at a particular point.
We saw in the Computational Stability article that a Fourier series technique could be used to determine a set of three stability conditions for the difference equation, Eq.2. Here we shall see what can be learned from looking at the truncation errors associated with the approximating equation, Eq. 2.
Truncation Error Evaluation
Assume that each term in Eq. 2 is a continuous and differentiable function of x and t. Then, for example, “uj+1,n would be u(xj+δx,tn) and can be expanded about the point (xj,tn) in a Taylor series in powers of δx. Carrying out the expansion in δx and δt for all the terms in Eq. 2 yields,
(3) $latex \displaystyle \frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}-\nu \frac{{{\partial }^{2}}u}{\partial {{x}^{2}}}=-\frac{1}{2}\delta t\frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}+O\left( \delta {{x}^{2}},\delta {{t}^{2}} \right).$
All second and higher order terms in δx and δt have been lumped into the order symbol O(δx2 ,δt2). This is a consistent approximation because it reduces to the original partial differential equation, Eq. 1, when δx and δt tend to zero.
Comparison of Fourier and Truncation Error Analysis
In the article Computational Stability a typical Fourier mode of the form
$latex \displaystyle P_{j}^{n}\propto {{r}^{n}}{{e}^{{ikxj}}}$ was substituted into the difference equation, Eq. 2, to obtain an equation for r,(4) $latex \displaystyle r=1-\left( \frac{ic\delta t}{\delta x} \right)\sin \left( k\delta x \right)-\left( \frac{2\nu \delta t}{\delta {{x}^{2}}} \right)\left[ 1-\cos \left( k\delta x \right) \right].$
Computational stability of the difference equation requires that the magnitude of r remain less than or equal to 1.0.
If we insert a Fourier mode of the form exp(i(kx+wt)) into the truncated Eq. 3, it will be seen that the result is the same as Eq. 4 with r=exp(iwδt) and then expanded in powers of wδt, plus the sine and cosine expanded in powers of kδx. This confirms that the two results are the same, as they should be to O(δx2,δt2) retained in Eq. 3.
However, the comparison also indicates that to keep the basic form of r in Eq. 4, with its real and imaginary parts, we must keep at least the first non-zero terms from the sine and cosine when they are expanded in powers of kδx. The first non-zero term in the imaginary contribution to r comes from sin(kδx) and is proportion to kδx, which corresponds to the first derivative with respect to x in Eq.3. The first non-zero term in the real part of r (other than 1) comes from cos(kδx) and is proportional to (kδx)2, which corresponds to the second derivative with respect to x in Eq. 3.
These observations lead to the rule-of-thumb that for the truncated equation to reproduce the lowest order real and imaginary parts of the amplification factor r, it is necessary to retain the lowest order even and odd derivatives with respect to each independent variable in the truncation error. In Eq. 3 there is only one first order term proportional to δt and it is a second derivative with respect to t. There are no first order terms proportional to δx.
Examining the Truncated Equation for Stability
Using the above rule-of-thumb, the truncated equation is,
(5) $latex \displaystyle \frac{\delta t}{2}\frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}+\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}-\nu \frac{{{\partial }^{2}}u}{\partial {{x}^{2}}}=0$
The first important thing to note is that this is not identical to the original partial differential equation, Eq. 1. The claim made here is that Eq. 5 is a better approximation of the finite-difference equation than Eq. 1 and because of this we can obtain information about the stability properties of the difference equation. This, in fact, is the case.
Recall that the difference equation propagated information into a region of influence bounded by lines whose slopes are dx/dt=±δx/δt. Similarly, the truncated Eq. 5 has a hyperbolic (i.e., wave) character because of the second space and second time derivatives, and the effective wave speeds are ±(2ν/δt)½. If the difference equation is to have any hope of approximating the truncated equation then its region of influence must at least encompass the region of influence of the truncated equation, which leads to the condition
(6) $latex \displaystyle \frac{2\nu }{\delta t}\le {{\left( \frac{\delta x}{\delta t} \right)}^{2}}$ or $latex \displaystyle \frac{2\nu \delta t}{\delta {{x}^{2}}}\le 1.$
Courant, Friedrichs and Lewy [2] used a similar region of influence condition, now called the Courant condition, which restricts the distance a wave travels in one time increment to less than one space increment. A violation of the Courant condition leads to an oscillating and exponentially growing instability. Condition Eq. 6 is precisely one of the stability conditions found from Fourier analysis in Computational Stability.
A similar Courant-type condition can be inferred from the two first order derivative terms (the advective terms) in the truncated Eq. 5, which propagate information with speed c,
(7) $latex \displaystyle \frac{c\delta t}{\delta x}\le 1.$
This stability condition, also identified in Computational Stability, likewise leads to an oscillating and growing instability when violated.
To uncover a third stability condition we must first rewrite the truncated equation by converting the δt term to have space instead of time derivatives, but in a way that still maintains the first order of the expansion. This is done by differentiating Eq. 3 by t and neglecting all first and higher order terms,
(8) $latex \displaystyle \frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}+c\frac{\partial }{\partial x}\frac{\partial u}{\partial t}-\nu \frac{{{\partial }^{2}}}{\partial {{x}^{2}}}\frac{\partial u}{\partial t}=O\left( \delta t \right)$
Next replace the first time derivative of u by t in this equation using Eq. 1 to obtain
(9) $latex \displaystyle \frac{{{\partial }^{2}}u}{\partial {{t}^{2}}}={{c}^{2}}\frac{{{\partial }^{2}}u}{\partial {{x}^{2}}}-2c\nu \frac{{{\partial }^{3}}u}{\partial {{x}^{3}}}+{{\nu }^{2}}\frac{{{\partial }^{4}}u}{\partial {{x}^{4}}}+O\left( \delta t \right)$
Finally, rewrite the truncated Eq. 5 using this result for the δt term
(10) $latex \displaystyle \frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=\left( \nu -\frac{{{c}^{2}}\delta t}{2} \right)\frac{{{\partial }^{2}}u}{\partial {{x}^{2}}}+c\nu \delta t\frac{{{\partial }^{3}}u}{\partial {{x}^{3}}}-\frac{{{\nu }^{2}}\delta t}{2}\frac{{{\partial }^{4}}u}{\partial {{x}^{4}}}.$
This result is identical to what would have been obtained by Taylor expanding the original finite-difference equation about the point x=jδx and t=(n+½)δt (and would probably have been easier).
According to our rule-of-thumb the last two terms on the right side proportional to δt can be dropped because they involve higher order derivatives than what is in the first δt term on the right side, which leaves,
(11) $latex \displaystyle \frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=\left( \nu -\frac{{{c}^{2}}\delta t}{2} \right)\frac{{{\partial }^{2}}u}{\partial {{x}^{2}}}.$
This is an alternative form for the truncated equation that retains only the lowest order (first) truncation errors and only those that contain the lowest even and odd derivatives with respect to each independent variable.
Equation 11 is nearly the same as the original Eq. 1, except for a modified diffusion coefficient. The significant thing here is that the diffusion coefficient can be negative. As long as the diffusion coefficient is positive solutions of Eq. 11 exhibit exponentially damped behavior, but with a negative coefficient solutions have an exponentially growing character, i.e., a computational instability! Thus, a further condition for computational stability is that the diffusion coefficient remains positive,
(12) $latex \displaystyle \frac{{{c}^{2}}\delta t}{2}\le \nu$
In this case the instability is a pure growing one without the oscillations in sign associated with the two earlier region-of-influence conditions. If instability is encountered, knowing whether it is exhibiting an oscillation in sign or not will identify it as either a region-of-influence violation or a negative diffusion coefficient. Having this knowledge makes it easier to find a remedy for the instability.
Application to Two-Dimensional Fluid Flow
A two-dimensional example (x,z) of water flowing under a laboratory scale sluice gate offers a test for examining a computational instability arising from non-linearity in the governing equations. The physical problem consists of water held behind a gate with an elevation of 0.9ft. Downstream (right) of the gate there is a water pool of depth 0.14 ft. Gravity is 32.2 ft/s2 in the negative z direction (down). At time t=0 the gate is raised up a distance of 0.125ft and water surges out into the pool. Figure 1 shows the resulting flow obtained with a Navier-Stokes solver [3] at t=0.35s. The solver used for this example has been optimized to automatically eliminate instabilities so none are apparent in this case, but it is possible to force the program to use non-optimum settings.
To demonstrate some unstable behavior we first examine a heuristic analysis performed on the vertical velocity equation used in the simulation. Focus is on the effective diffusion coefficients for the z direction velocity w, while all other truncation errors are ignored,
(13) $latex \displaystyle \frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+w\frac{\partial w}{\partial z}+\frac{\partial }{\partial z}\left( \frac{p}{\rho } \right)+g=\left( \nu +\frac{\alpha u\delta x}{2}-\frac{{{u}^{3}}\delta t}{2}-\frac{\delta {{x}^{2}}}{4}\frac{\partial u}{\partial x} \right)\frac{{{\partial }^{2}}w}{\partial {{x}^{2}}}+\left( \nu +\frac{\alpha w\delta z}{2}-\frac{{{w}^{2}}\delta t}{2}-\frac{\delta {{z}^{2}}}{2}\frac{\partial w}{\partial z} \right)\frac{{{\partial }^{2}}w}{\partial {{z}^{2}}}$
The diffusion of w in the x and z directions are expressed by the two terms on the right side of Eq. 13, where ν is the fluid viscosity and α is a parameter that modifies the numerical approximation of the term describing the u advection of w, i.e., the second term on the left side of the above equation. When α=0 the finite-difference advection approximation is said to be centered about the location of w, but when α=1 an upstream or “donor cell” approximation is used.
The first thing to notice is that if ν=0 and a centered difference approximation is also used (α=0) then the lowest order term in the two effective viscosity coefficients are proportional to δt and are negative. This clearly leads to unstable behavior, and is a well known property of the central difference approximation. Adding enough viscosity to keep the diffusion coefficient positive is also an established procedure to gain stability, but at the possible cost of introducing too much diffusion. The upstream difference option, α=1, is a reasonable compromise; provided the condition wδt<δx is maintained, the diffusion coefficients are positive (provided the δx2 and δz2 terms are small) and the simulation will be stable.
If the δx2 and δz2 terms in the diffusion coefficients are not small there is a possibility of unstable behavior. To demonstrate this we set the viscosity to zero and reduce the amount of upstream differencing by setting α=0.05. To keep the negative δt term less than the a term a very small time step δt=0.00025 is used. With these settings the resulting simulation is shown in Fig. 2. An instability in the z velocity has developed just upstream of the sluice gate, which is shown close up in Fig. 3 (where color indicates the z velocity magnitude).
This instability is a result of a negative x-direction diffusion coefficient, which is coming from the δx2 term. A negative value results from the fact that the flow upstream of the gate is compressing in the z direction, but expanding in the x direction, which means that the x derivative of u in the δx2 term is positive in this region resulting in a net negative diffusion coefficient.
A check on this conclusion can be made by adding in a little viscosity ν=0.0093 to compensate for the negative δx2 term. Figure 4 shows that this change does, indeed, stabilize the flow.
This example demonstrates that truncation error terms arising from non-linear terms in the original equation influence the computational stability of the difference equation. This type of instability cannot be found by a von Neumann type Fourier analysis. Perhaps most important of all is that when troublesome truncation errors are found to exist this knowledge can be used to alter the finite difference equations to eliminate those errors.
Summary
To summarize, it has been shown that all the stability conditions associated with a linear finite-difference equation, Eq.2, can be identified using a heuristic truncation error approach. This approach not only identifies the instabilities, it also indicates what can be done to eliminate them. For instance, for a region-of-influence violation only a reduction in the time-step increment will solve the problem, but if there is a negative diffusion coefficient then adding more diffusion to compensate for the errors is one way to regain stability. Knowing the origin of a negative diffusion error may also suggest how the original finite-difference equation might be modified to avoid this problem.
The most significant aspect of the heuristic approach is that it is not limited to linear equations with constant coefficients, as was shown in connection with the example of flow under a sluice gate. No special assumptions were necessary to form the approximating truncated equation. The goal was simply to reverse the procedure of writing a difference equation to approximate a partial differential equation, and instead to write a partial differential equation that approximates the difference equation. A simple rule-of-thumb was described for constructing the truncated equation. This approximating equation was then used to check for region-of-influence violations and for possible negative diffusion coefficients both features that lead to unstable solutions.
Several additional examples involving compressible and incompressible fluid dynamics simulations can be found in the original heuristic stability paper [1], which further show how the heuristic approach can be applied to real, practical, non-linear problems.
