Box 6.1 Stiff Systems of Equations
A system of first-order differential equations in time
(I)
is said to be stiff if the timescales for change of the dependent variables ψi range over many orders of magnitude. We saw in Section 4.4 that the characteristic timescales of the system are defined by the inverses of the eigenvalues λi of the Jacobian matrix J = ∂s/∂ψ where s = (s1, … sK)T and ψ = (ψ1, … ψK)T. Stability of the solution requires that the real component of the eigenvalues be negative and this is always met in relevant mechanisms for atmospheric chemistry (with only transient exceptions). Thus
Re(λi) < 0 for i = 1, 2, . . . K
(II)
We define the stiffness ratio R as
(III)
A simple example of a stiff system is given by (Press et al., 2007):
(IV)
with the initial conditions u(0) = 1 and v(0) = 0. The two eigenvalues of this system are –1 and –1000, corresponding to characteristic timescales of 1 and 0.001, and a stiffness ratio of 103. The analytical solution is:
(V)
Box 6.1 Figure 1 shows the exact solution for u for t ∈ [0, 0.02] as well as numerical solutions by fully explicit and implicit algorithms with different time steps. u increases from 1 to 2 over the first characteristic timescale of 0.001, and then declines on the much longer timescale of 1 (not shown in the figure). The implicit algorithm is stable for all time steps and provides accurate solution of the asymptotic behavior for the first characteristic timescale, though not of the transient behavior. The explicit algorithm is stable when ∆t = 0.001 but incurs oscillatory behavior at longer time steps that becomes undamped for ∆t = 0.002. Even though the e–1000t term becomes rapidly negligible as t increases, the explicit solution still requires a time step of less than 0.002 throughout the integration. This is burdensome if we are interested in computing the solution over long times (t ≫ 1) in order to capture the second characteristic timescale of the system.
Box 6.1 Figure 1 Solution of stiff system (IV) for the range 0 < t < 0.02 and three values of the time step ∆t. The exact analytical solution for u is compared to numerical solutions by fully implicit and explicit methods and for different time steps.
An alternate approach is to use an implicit solver, where s in (6.1) is estimated on the basis of the unknown concentrations at time to + Δt. The system of coupled ODEs then becomes a system of K coupled algebraic equations to solve for ψ(to + Δt). Implicit methods have far less severe restrictions on size of time step to remain stable. However, they require computationally expensive constructions to obtain the solution of the system of equations. High-order implicit solvers such as Rosenbrock and Gear are often used as standards of accuracy in 3-D models. Other algorithms provide a compromise between accuracy and computational performance.
The choice of a particular numerical algorithm for the chemistry operator depends on several considerations including stability, positivity, accuracy, mass conservation, and computational efficiency of the method (Zhang et al., 2011). Positivity of the solution is essential for stability as otherwise the kinetics equations immediately diverge. This condition can severely restrict the size of the time step. Mass conservation is essential if quantitative tracking of chemical budgets is needed. Some slack in accuracy is often considered acceptable because the kinetics equations have stable solutions (Section 4.4) so that inaccuracies will dampen rather than grow. Computational efficiency may be of no concern for a box model but critically important for a 3-D model.
6.2 General Considerations
6.2.1 Fully Explicit Equation
The simplest method to solve (6.1) is the single step forward Euler or fully explicit algorithm:
ψn + 1 = ψn + Δt s(tn, ψn)
(6.3)
where n is a time index, tn = t0 + n ∆t (n = 0, 1, 2, …) is the discretized time level, ∆t is the integration time step, and ψn and ψn+1 are approximations to the solution ψ(t) at time levels tn and tn+1, respectively. The source term s(tn, ψn), also noted sn, is evaluated at time tn for the known approximation ψn. This method, which is first-order accurate, is called fully explicit because the unknown ψn+1 is represented strictly as a function of the known quantities at previous time tn. It is mass-conserving but positivity is not guaranteed. Equation (6.3) is a single step algorithm because only time levels tn and tn+1 are involved. Higher-order multi-step fully explicit schemes that express ψn+1 as a function of the solution at previous time levels tn, tn–1, tn–2, etc. are described in Section 6.2.3.
The explicit algorithm is appealingly simple but suffers from severe stability restrictions, as can be illustrated with a trivial example. Consider a single chemical species subject only to a linear loss. The corresponding differential equation is
(6.4)
with loss coefficient ℓ assumed to be constant. Applying the forward Euler algorithm, we have
ψn + 1 = ψn(1 − ℓΔt)
(6.5)
which is an approximation to the exact analytic solution
ψ(tn + 1) = ψ(tn) exp [−ℓΔt]
(6.6)
Stability requires that ⃒ψn+1/ψn⃒ < 1 as the integration proceeds. Equation (6.5) meets this stability criterion if
(6.7)
Positivity of the solution requires the more stringent criterion
(6.8)
Thus the time step for a system of several chemical species must be smaller than the chemical lifetime of the fastest reacting species.
6.2.2 Fully Implicit Equation
The constraint on the time step associated with fully explicit methods can be overcome by using the backward Euler or fully implicit algorithm. In this case, the solution to (6.1) is approximated by
ψn + 1 = ψn + Δt s(tn + 1, ψn + 1)
(6.9)
in which the source term s is evaluated at time tn+1 and is therefore expressed as a function of the unknown quantity ψn+1. The fully implicit algorithm ensures the positivity of the solution and also conserves mass. The resulting system of algebraic equations in ψn+1 requires numerical solution except in trivial cases.
In the simple linear example (6.4), the backward Euler scheme leads to
(6.10)
We have ⃒ψn+1/ψn⃒ < 1 for any positive value of ∆t, thus the numerical scheme is unconditionally stable. The ratio tends to zero for large values of ∆t, mirroring the analytic solution. The algorithm is first-order accurate; the error in the solution can be estimated by comparing the analytic solution exp [–ℓ ∆t] with the approximation 1/(1 + ℓ ∆t).
Table 6.1 shows the solution of (6.4) for t = 2, when ℓ is chosen to be 1 and the initial value is ψ(0) = 1. The exact (analytic) solution is compared to the approximate solution derived with the fully explicit and fully implicit algorithms for different values of the time step ∆t. As expected, for both algorithms, the accuracy of the solution decreases as the time step increases. Although the implicit algorithm has the advantage of stability, that does not make it any more accurate.
Table 6.1 Solution of equation dψ/dt = –ψ at time t = 2 with ψ(0) = 1
Solution at time t = 2
Time step Exact Fully explicit Fully implicit
0.0001 0.1353 0.1353 0.1353
0.001 0.1353 0.1352 0.1355
0.01 0.1353 0.1340 0.1367
0.1 0.1353 0.1216 0.1486
1.0 0.1353 0.0000 0.2500
2.0 0.1353 –1.0000 0.3333
A simple approach to solve implicit equation (6.9) is to linearize the source term s around the solution ψ at time tn:
s(tn + 1, ψn + 1) = s(tn, ψn) + J(ψn + 1 − ψn)
(6.11)
where J = ∂s/∂ψ is the Jacobian matrix of partial derivatives estimated for ψ = ψn with elements Ji,j = ∂si/∂ψj. Using (6.11) in (6.9) we obtain
ψn + 1 = ψn + Δt [I − JΔt]−1s(tn, ψn)
(6.12)
in which I is the identity matrix. Applying this approach involves solving a matrix system of equations. This li
nearization method, when applied to an implicit equation, is called the semi-implicit Euler method. It is usually but not unconditionally stable. Other methods to solve implicit equations are discussed in Section 6.4.
6.2.3 Improving Accuracy
Both the Euler forward and backward methods are asymmetrical since the time derivatives are evaluated in one case at the beginning of the time interval and in the other case at the end of the interval. They are therefore only first-order accurate in ∆t. Figure 6.1 illustrates the difference between the forward and backward methods, highlighting the errors incurred in the first-order approximation.
Figure 6.1 Determination of ψ(t) at time level tn+1 from its known value at time tn. The tangent A–B1 is proportional to the source term s at time tn+1, while the tangent A–B2 is proportional to the source term at time tn. Points B1 and B2 represent approximations to the true solution B obtained by the implicit and explicit Euler algorithms, respectively. The respective errors are defined by the distances B–B1 and B–B2.
The accuracy of the solution ψn+1 can be improved by making the solver more symmetric relative to time levels tn and tn+1. This can be done by taking the average of s between time levels tn and tn+1, which is equivalent to adopting a time-centered derivative:
(6.13)
Equation (6.13) defines the Crank–Nicholson scheme. This semi-implicit algorithm is second-order accurate. Numerical solution is required as in the backward Euler fully implicit scheme. The solution is not guaranteed to be positive.
It is also possible to increase accuracy in the framework of an explicit algorithm by using predictor-corrector methods. In such a method, the prediction step derives a first estimate of the solution (un+1) at time tn+1 from the forward Euler equation:
un + 1 = ψn + k1Δt
(6.14)
where the slope k1 = s(tn, ψn) is calculated at time level tn. The solution is improved by applying a correction step in which the slope is replaced by the average of k1 at time tn and an estimate k2 = s(tn+1, un+1) at time tn+1:
(6.15)
It can be shown that if the third derivative of the solution is a continuous function, this improved Euler method is a second-order scheme. The choice of the time step ∆t remains constrained by the stability criteria of explicit methods.
In the midpoint method, the solution ψn+1 is derived from a Euler formula in which s is estimated at an intermediate time level tn+1/2 = tn + ∆t/2. In this algorithm, the first step is to derive an estimate un+1/2 of the solution at midpoint of interval ∆t
(6.16)
with again k1 = s(tn, ψn). In the second step, the solution is computed using the entire time interval
ψn + 1 = ψn + k2Δt
(6.17)
where k2 = s(tn+1/2, un+1/2) is an estimate of the source term at the midpoint between time levels tn and tn+1. Due to its symmetrical nature, the midpoint method is second-order accurate.
The accuracy of the solution can also be improved by applying an s-stage Runge–Kutta method defined by
(6.18)
where
(6.19)
with bi(i = 1, … s) and ai,j(i, j = 1, … s) chosen to meet desired accuracy and stability conditions, and ci = Σj ai,j. If all coefficients ai,j ≠ 0, the method is fully implicit and numerically highly stable. Most applications, however, are based on the explicit Runge–Kutta method in which coefficients ai,j = 0 for j ≥ i. The method is less robust than the implicit version but is easier to apply. In the case of the classic explicit fourth-order Runge–Kutta method, for example, the solution at time tn+1 is provided by
(6.20)
where
k1 = s(tn, ψn) represents the source term at time level tn,
k2 = s(tn + ∆t/2, ψn + k1 ∆t/2) denotes a first estimate of the source term at the midpoint of the interval [tn, tn+1],
k3 = s(tn + ∆t/2, ψn + k2 ∆t/2) represents an improved estimate of the source term at the midpoint, and
k4 = s(tn+1, ψn + k3 ∆t) estimates the source term at time level tn+1, using the value of k3 calculated at the midpoint.
Explicit Runge–Kutta methods are more stable than the forward Euler algorithm. They are usually implemented with an adaptive step size procedure to meet a user-required error tolerance. As in other explicit methods, the time step is constrained by the shortest lifetime in the system and this can make implementation for atmospheric problems impractical. Implicit Runge–Kutta methods (Hairer et al., 2002) are characterized by high stability, but the resulting system of equations is difficult to solve. The diagonally implicit Runge–Kutta method (DIRK), in which ai,j = 0 for j > i, but the diagonal elements aii ≠ 0, is simpler to implement than the fully implicit case. The RADAU5 solver, introduced by Hairer et al. (1993) and implemented in some chemical models, is a one-step implicit Runge–Kutta method that is fifth-order accurate.
Linear multi-step algorithms, which retain the information from several previous time steps (tn, tn–1, tn–2, etc.), also provide solutions with higher order of accuracy. They can be expressed either by an explicit expression
(6.21)
or by an implicit expression
(6.22)
where bj are constant coefficients and s corresponds to the order of the method. When combined, these two relations represent a multi-step predictor-corrector scheme. In this case, function s(tn+1, ψn+1) appearing in the correction step (6.22) is estimated by using the values of ψn+1 derived by the prediction step (6.21). An example is the Adams–Bashforth–Moulton scheme. For a value s = 3, for example, the prediction step is
(6.23)
and the correction step is
(6.24)
Even though (6.22) is an implicit expression, the introduction of a predictor-corrector approach transforms the scheme into a fully explicit scheme with the associated stability requirements.
Accuracy can also be improved by applying the extrapolation method introduced by Lewis Fry Richardson. This method allows the construction of high-order solutions by applying the same algorithm with decreasing time steps. It is based on asymptotic expansion of the truncation error in h-powers, where h is the time step. The approximate solution ψn+1(h) at time tn+1 derived with some numerical algorithm Lh and time step h can be expressed as
(6.25)
where ψ(tn+1) represents the true solution at time tn+1. Index m represents the order of the scheme, while Em(ψ)hm and O(hm + 1) are the errors associated with algorithm Lh of order hm and hm+1, respectively. When the same algorithm is applied with smaller time steps h/k (k = 1,2,3, …), we write similarly
(6.26)
Combining (6.25) and (6.26) yields Richardson’s recurrence formula that provides a higher-order approximation
(6.27)
If, for example, k = 2 and if ψn+1(h) and ψn+1(h/2) are the numerical solutions obtained by a first-order algorithm (m = 1) with time steps h and h/2, respectively, the accuracy of the solution is improved by using
ψn + 1 = 2 ψn + 1(h/2) − ψn + 1(h)
(6.28)
Examples of extrapolation methods are given in Sections 6.3.1 and 6.3.3.
6.2.4 Explicit Versus Implicit Solvers
The comparison between fully explicit and fully implicit methods highlights the advantages and disadvantages of both approaches (Sandu et al., 1997b). Fully explicit equations are usually simple to solve, but stability and positivity considerations may constrain the integration time steps to prohibitively small values. Fully implicit methods are unconditionally stable and positive, so that the time step can be large, limited by accuracy requirements. However, they require solution of a system of algebraic equations at each time step, involving in general the construction and inversion of a Jacobian matrix. This can be computationally costly. Methods have been developed to reduce the stiffness of chemical systems by separating species between short-lived and long-lived and solving for each group separately, with an implicit method used for the short-lived subset only (Gong and Cho, 1993). However, separation is often difficult because at
mospheric chemical mechanisms typically involve a continuum of lifetimes and the lifetimes vary with the local conditions. Adaptive separation can be done locally within a model simulation by calculating species lifetimes before applying the chemistry operator (Santillana et al., 2010), but this involves substantial computational overhead.
6.3 Explicit Solvers
We first present several explicit algorithms that use various methods to relax the requirement for short time steps while keeping the computational advantage of the explicit solution. We rewrite the system of equations to separate production and loss terms as
Modeling of Atmospheric Chemistry Page 33