Numerical Approximations for Partial Differential Equations
Table of Contents
Mathematical models for continuum dynamic phenomena involve one or more partial differential equations. Constructing numerical approximation for these equations that can be used to obtain approximate solutions using computers requires some care. For poorly selected approximations, the numerical results produced can be at best highly inaccurate or at worst unstable and completely unusable. To gain some understanding of how numerical approximations work, it is useful to study some simple examples whose details are easy to work out.
For this purpose, we will study the simple example of the first order transport equation of a scalar S with constant advection speed c in one dimension x,
$latex \displaystyle \frac{\partial S}{\partial t}+c\frac{\partial S}{\partial x}=0$ (1)
Deriving a Numerical Difference Approximation
A numerical, finite difference, solution is to be constructed on a grid of points on the x axis that are spaced a distance δx apart. Location xj would be at the jth δx. The notation Sjn is the value of S at x location jδx and time t=nδt.
Each derivative in Eq. 1 can be approximated using a Taylor series expansion in powers of δt and δx. For example,
$latex \displaystyle S_{j}^{n+1}=S+\delta t\frac{\partial S}{\partial t}+\frac{\delta {{t}^{2}}}{2}\frac{{{\partial }^{2}}S}{\partial {{t}^{2}}}+…$
and
$latex \displaystyle S_{j+1}^{n}=S+\delta x\frac{\partial S}{\partial x}+\frac{\delta {{x}^{2}}}{2}\frac{{{\partial }^{2}}S}{\partial {{x}^{2}}}+…$ (2)
Rearranging in terms of the first derivative, a forward difference for ∂S/∂t and central difference for ∂S/∂x can be derived as follows:
$latex \displaystyle \frac{\partial S}{\partial t}=\frac{S_{j}^{n+1}-S_{j}^{n}}{\delta t}+O(\delta t)$
and
$latex \displaystyle \frac{\partial S}{\partial x}=\frac{S_{j+1}^{n}-S_{j-1}^{n}}{2\delta x}+O(\delta {{x}^{2}})$ (3)
Substituting Eq. 3 back into Eq. 1, we arrive at a simple numerical approximation of a forward difference for ∂S/∂t and central difference for ∂S/∂x:
$latex \displaystyle S_{j}^{n+1}-S_{j}^{n}=-\left( \frac{c\delta t}{2\delta x} \right)\left( S_{j+1}^{n}-S_{j-1}^{n} \right)+O(\delta {{t}^{2}},\delta {{x}^{2}})$ (4)
It is evident that we cannot use an infinite number of terms from the Taylor series to form our numerical approximations, and therefore must neglect these higher order terms. The terms O(δt2) and O(δx2) represent the leading terms truncated from each approximation and may also be referred to as the leading error terms. The exponent of the leading error term defines the formal order of accuracy and gives a measure of how rapidly the error will improve with additional refinement of the discrete interval.
Numerical Stability
A critical aspect of numerical solutions is that they can become unstable by growing exponentially. When this behavior occurs, oscillating features will be observed that bear no relation to the solution of the original differential equation. If we want to understand whether a numerical approximation is stable, a customary method is to use a Fourier series approach first devised by mathematician John von Neumann1.
One assumes the solution of the difference equation, Eq. 4, has the form for the kth Fourier mode with amplitude r
$latex \displaystyle S_{j}^{n}=\tilde{S}{{r}^{n}}{{e}^{ikj\delta x}}$ (5)
Substituting this into Eq. 4 (where i equals the square root of -1) gives
$latex \displaystyle r=1-\left( \frac{c\delta t}{\delta x} \right)i\sin \left( k\delta x \right)$ (6)
Using the following relation
$latex \displaystyle a+ib=\sqrt{{{a}^{2}}+{{b}^{2}}}{{e}^{i{{\tan }^{-1}}\left( b/a \right)}}$ (7)
and assuming cδt/δx is small the kth mode solution is
$latex \displaystyle S_{j}^{n}=\tilde{S}{{A}^{n}}{{e}^{ik\left( x-{c}^{*}t \right)}}$ (8)
where
$latex \displaystyle {{A}^{2}}=1+{{\left( \frac{c\delta t}{\delta x}\sin \left[ k\delta x \right] \right)}^{2}}$ and $latex \displaystyle c*=c\frac{\sin \left( k\delta x \right)}{k\delta x}$
What this reveals is that the amplitude A is greater than 1.0 which means the amplitude of the mode is growing by the power of n, so the approximation is unconditionally unstable and will not produce any meaningful result. A secondary finding is that the speed of the Fourier mode, c*, depends on the wave number k and this means that different wave components of the solution will not maintain their summation causing oscillations to appear in the solution, a phenomenon called dispersion.
Another useful analysis tool is to look at the truncation errors associated with the approximation of Eq. 2. For this purpose, each term in the equation is expanded in a Taylor series in powers of δt and δx, see Eq. 2. When applied to Eq. 4 the result to first order in δt and δx, is
$latex \displaystyle \frac{\partial S}{\partial t}+c\frac{\partial S}{\partial x}=-\frac{\delta t}{2}\frac{{{\partial }^{2}}S}{\partial {{t}^{2}}}+higherorderterms$ (9)
This appears to be second order in time because of the second order time derivative, but that would require two initial conditions when we know that only one initial condition is needed, so this is not a proper or consistent approximation of the difference equation. This can be resolved by differentiating the zeroth order terms by t and replacing the second time derivative by a second space derivative to arrive at
$latex \displaystyle \frac{\partial S}{\partial t}+c\frac{\partial S}{\partial x}=-\frac{{{c}^{2}}\delta t}{2}\frac{{{\partial }^{2}}S}{\partial {{x}^{2}}}+higherorderterms$ (10)
Ignoring the higher order terms shows that the truncation error term on the right side of the equation appears to be introducing diffusion, but the coefficient is negative. A negative diffusion implies a growth of S variations and not decay. This is another way to find that the difference approximation of Eq. 4 is unconditionally unstable. Correspondingly, a positive diffusion would be stable, but would introduce some damping in the values of S.
Numerical Error
Before moving on it is worthwhile to note a couple of things about Eq. 10. For instance, if the space and time increments go to zero, the equation reduces to the original partial differential equation. Because of this, the difference approximation is said to be consistent. With nonzero increments δx and δt, however, the approximation contains two new parameters that are not in the original equation. Dimensional analysis would suggest that the approximate equation should have a greater range of solutions than the original. The low order truncation errors do add some diffusion (positive and negative) that is not in the original formulation and this diffusion can be seen in the computational results. For this reason, Eq. 10 is a more accurate partial differential equation representing the difference equation and is sometimes referred to as the modified equation 2.
Other Difference Approximations
Having discussed two methods to check the stability of a finite difference approximation, the next step is to consider alternative approximations and to learn something more about the accuracy and stability of numerical approximations. In addition to the above analysis, computational results will be shown for a simple test problem. For this test, the advection of the scalar variable S will be computed. The initial value of S is a unit step function consisting of S=1 in the left half of the computational grid and 0 in the remaining half, see Fig. 1. The advection velocity will be a constant of 1.0 from left to right. Since S is constant on either side of the step the only interesting region to observe changes in S is around the step.
In the following examples several different approximations for advection will be used to illustrate how their stability and dispersion properties can affect the results they produce. A short computer program has been used to solve this numerical test problem. This program is listed in the Appendix in Small Basic and Python. Readers are encouraged to try one of these programs and select input for different numerical approximations to see the results they produce.
Example 1: Central Difference Advection
The sample difference equation in Eq. 4 is referred to a centered difference approximation because the advection term derivative is approximated by the difference of S values centered about the Sj
$latex \displaystyle S_{j}^{n+1}-S_{j}^{n}=-\left( \frac{c\delta t}{2\delta x} \right)\left( S_{j+1}^{n}-S_{j-1}^{n} \right)$ (11)
It has already been shown that the amplitude A of the Fourier components making up S is greater than unity,
$latex \displaystyle {{A}^{2}}=1+{{\left( \frac{c\delta t}{\delta x}\sin \left[ k\delta x \right] \right)}^{2}}$ (12)
which means that as the time step n increases, the amplitude increases exponentially. It was also noted that the phase of the Fourier modes depends on the wave number k, and this means there will be dispersion causing the sum of modes to not be able to combine into fixed values (i.e., spatial oscillations will appear in S).
The truncation equation for Eq. 10 was shown to have a first order in time truncation error,
$latex \displaystyle \frac{\partial S}{\partial t}+c\frac{\partial S}{\partial x}=-\frac{{{c}^{2}}\delta t}{2}\frac{{{\partial }^{2}}S}{\partial {{x}^{2}}}+O\left( \delta {{t}^{2}},\delta {{x}^{2}} \right)$ (13)
As previously noted, the error on the right-hand side looks like a diffusion but really has a negative diffusion coefficient, which indicates that a solution of this equation will show exponentially growing instabilities.
A numerical solution of Eq. 10 after 25 time steps, where cδt/δx=0.5, is shown in Fig. 1 and confirms the expected behavior that was predicted by both the Fourier analysis and truncation error analysis of the difference equation.
Before progressing to other examples, the nondimensional combination of parameters cδt/δx appears in Eq. 11. This quantity is called the Courant Number, which is the ratio of the distance that advection moves S values in one time step δt divided by the width of the mesh cell δx.
The Courant Number plays a prominent role in problems involving advection. From a physical perspective it should be apparent that moving S values by more than one mesh interval δx cannot be accurate because the numerical approximation in Eq. 10 only relates changes in Sj caused by its immediate neighbors Sj-1 and Sl+1. Information about S values further away than one δx from Sj have no influence on Sj. For the remainder of this document, it will be assumed that the Courant Number is always less than or equal to 1.0. Thus, there will be no issues with a possible loss of accuracy because of an insufficient region of influence in the approximations.
Example 2: Implicit Central Difference Advection
In this approximation the only change from the former case is the advanced time level n+1 in the advection term. The presence of n+1 values on the right-hand side of the equation, which are initially unknown makes the equation implicit. In general, implicit equations are much more difficult to solve, but in this one-dimensional case there is a simple analytic solution.
$latex \displaystyle S_{j}^{n+1}-S_{j}^{n}=-\left( \frac{c\delta t}{2\delta x} \right)\left( S_{j+1}^{n+1}-S_{j-1}^{n+1} \right)$ (14)
A Fourier stability analysis yields the mode amplitude,
$latex \displaystyle {{A}^{2}}=\frac{1}{1+{{\left( \frac{c\delta t}{\delta x}\sin \left[ k\delta x \right] \right)}^{2}}}$ (15)
Where it is clear that A is always less than 1.0, so the numerical scheme has no growing modes and is therefore stable. Furthermore, because the wave number k appears in the imaginary portion of the amplitude there will be dispersion.
A truncation analysis of Eq. 13 results in,
$latex \displaystyle \frac{\partial S}{\partial t}+c\frac{\partial S}{\partial x}=\frac{{{c}^{2}}\delta t}{2}\frac{{{\partial }^{2}}S}{\partial {{x}^{2}}}+O\left( \delta {{t}^{2}},\delta {{x}^{2}} \right)$ (16)
Here, as in Eq. 13, there is a diffusion-like term on the righthand side of the equation, but now the diffusion coefficient is positive, meaning we should expect stability, but also some damping or dissipation to be evident.
A plot of the numerical solution produced by this approximation after 40 time steps is shown in Fig. 2. In this case, the Courant Number is 0.5.
Because the mode amplitude is always less than 1.0, the consequences of this are a significant damping or smoothing of the initial step in S that is moving to the right (this was expected from the truncation error). Some dispersion effects, arising from the initial step in S, can also be seen traveling to the left.
Example 3: Upstream or Donor Differencing for Advection
A popular difference method is to compute advective fluxes using upstream values. For example, using the upstream jth value of S instead of j+½ and the upstream j-1 value instead of j-½. Since the advection is from left to right the flux across a cell boundary would seem to be more physical if weighted upstream to the left.
$latex \displaystyle S_{j}^{n+1}-S_{j}^{n}=-\left( \frac{c\delta t}{\delta x} \right)\left( S_{j}^{n}-S_{j-1}^{n} \right)$ (17)
A Fourier stability analysis produces the mode amplitude,
$latex \displaystyle {{A}^{2}}=1-4{{N}_{c}}(1-{{N}_{c}}){{\sin }^{2}}\left( \frac{k\delta x}{2} \right)$ (18)
Here Nc is the Courant Number cδt/δx. Amplitude A remains less or equal to 1.0 for all modes and all Courant Numbers (less than or equal to 1.0). The amount of damping depends on the value of Nc. There is also a k dependent imaginary portion that indicates dispersion will be present as well.
A truncation error analysis produces,
$latex \displaystyle \frac{\partial S}{\partial t}+c\frac{\partial S}{\partial x}=\frac{c\delta x}{2}(1-{{N}_{c}})\frac{{{\partial }^{2}}S}{\partial {{x}^{2}}}+O\left( \delta {{t}^{2}},\delta {{x}^{2}} \right)$ (19)
Again there is a diffusion-like truncation error that remains positive as long as the Courant Number is not greater than 1.0. This supports the conclusion of the Fourier analysis that the numerical equation will exhibit some damping. However, in this case, the damping can be somewhat controlled by the value of the Courant Number such that a larger Courant Number, but not exceeding 1.0, will lead to less damping.
Numerical solutions of Eq. 17 are shown in Fig. 3 for three values of Courant Number.
As expected from the truncation error, when the Courant Number is 1.0 there is no damping of the step profile of S (left plot), while a Courant Number of 1.05 immediately leads to instability (middle plot). The right plot for a Courant Number of 0.5 is stable but exhibits a considerable amount of damping or smoothing of the step profile.
It is important to note that the effective diffusion coefficient in Eq. 18 is the difference between two truncation errors. The first term is cδx/2 and the second term is –c2δt/2. Each term is a first order truncation error, but the second one is negative. This is a case where there is competition between space and time errors and in particular where the δx term can compensate for the δt term, provided the Courant Number remains less than or equal to 1.0, to produce a stable difference approximation.
Example 4: Crank-Nicholson Differencing for Advection
A popular differencing method for higher order in time as well as space is to use centered advection but to split it half at the old time level n and half at the new time level n+1. Splitting the time level for advection ensures the time step truncation error will be second order.
$latex \displaystyle S_{j}^{n+1}-S_{j}^{n}=-\left( \frac{c\delta t}{4\delta x} \right)\left( S_{j+1}^{n+1}-S_{j-1}^{n+}+S_{j+1}^{n}-S_{j-1}^{n} \right)$ (20)
A Fourier stability analysis yields the mode amplification A to be,
$latex \displaystyle {{A}^{2}}=\frac{(1+N_{c}^{2})}{(1+N_{c}^{2})}=1$ (21)
This result implies there is no dissipation or damping, but dispersion is still expected. The method is stable for all Courant Numbers less than or equal to 1.0.
Figure 4 shows the numerical results obtained using this second order approximation for advection after 40 time steps. Dispersion errors are clearly evident and are larger than seen in previous approximations because there is no damping present.
Truncation errors are not displayed because they are zero through first order. Only second order errors are present. This is consistent with the Fourier stability analysis. That is, there is no stable or unstable first order diffusion-like truncation errors, so no growth or damping is seen in the numerical results. There are significant oscillations caused by dispersion effects, and without some damping, these oscillations pretty much spoil the computed results. A lesson to learn from this example is that going to a higher order approximation does not always produce better results.
Lax-Wendroff Approximation for Advection
The Lax-Wendroff numerical approximation is second order in both space and time. However, second order in time is arranged by a different method than used in the previous Crank-Nicholson method. The idea is to replace the second order in time truncation error by equivalent second order space derivatives, as computed from the original partial differential equation,
$latex \displaystyle S_{j}^{n+1}-S_{j}^{n}=-\left( \frac{c\delta t}{2\delta x} \right)\left( S_{j+1}^{n}-S_{j-1}^{n} \right)+\frac{1}{2}{{\left( \frac{c\delta t}{\delta x} \right)}^{2}}\left( S_{j+1}^{n}-2S_{j}^{n}+S_{j-1}^{n} \right)$ (22)
The Fourier mode amplitude for this numerical approximation is
$latex \displaystyle {{A}^{2}}=1-4N_{c}^{2}\left( 1-N_{c}^{2} \right){{\sin }^{4}}\left( \frac{k\delta x}{2} \right)$ (23)
This amplitude remains less than 1.0 (is stable) for Nc less than or equal to 1.0. The approximation is also dispersive. When Nc is close to 1.0 there is little or no damping because A~1.0. Numerical solutions verify this observation as seen in Fig. 5. When Nc=0.98, the computed results after 20 time steps (left plot) shows little smoothing of the initial step in S. At Nc=0.1, the computed results after 200 time steps (right plot) is beginning to show the dispersion seen in the Crank-Nicholson results, but with some damping of the oscillations. The Crank-Nicholson method has no damping so it exhibits somewhat larger oscillations.
Concluding Comments
The intent of this article is to show how the computational results of several finite difference approximations, while appearing different in many respects, can nevertheless be understood through simple analysis using both the Fourier and truncation error methods. Understanding these results in the simple test case used here can provide guidance in evaluating the results obtained in more complicated cases involving multiple equations and additional models for physical phenomena.
References
- See the article on Computational Stability in CFD-101
- For more about truncation error analysis see the Heuristic Analysis article in CFD-101.
Generating the Plots in this Document
The following program was written in both Microsoft® Small Basic and Python, which you can download below.
