3.5 Numerical solution of the split sub-problems
3.5.4 Model for a stiff problem: reaction-diffusion equation
When the operator splitting theory is applied to the numerical solution of partial differen-tial equations, it is usually assumed that in the first step we semi-discretize (in the space variable) the problem in order to get a Cauchy problem. Then we can use different split-tings. Consequently, for this case, the operators that appear in the split tasks depend on the space discretization parameter. Typically they result in a stiff problem, which, due to the stability, has special features and requires the use of special numerical methods [30].
Therefore it is worth considering the different operator splitting and numerical methods for such problems separately.
Consider the diffusion-reaction equations with linear reaction [73]:
∂u
∂t =D1∂2u
∂x2 −k1u+k2v+s1(x)
∂v
∂t =D2∂2v
∂x2 +k1u−k2v+s2(x),
(3.5.25)
for the unknown concentration functionsu(x, t) and v(x, t), where 0< x <1 and 0< t≤ T = 12, and the initial and boundary conditions are defined as follows:
u(x,0) = 1 + sin(0.5πx), v(x,0) = (k1/k2)u(x,0),
u(0, t) = 1, v(0, t) = k1/k2,
∂u
∂x(1, t) = ∂v
∂x(1, t) = 0.
(3.5.26)
We used the following parameter values:
• diffusion coefficients: D1 = 0.1; D2 = 0,
• reaction rates: k1 = 1; k2 = 104,
• source terms: s1(x)≡1; s2(x)≡0.
After second-order finite-difference space discretization, a two-stage diagonally implicit Runge–Kutta (DIRK) [74] method was used for time discretization. Table 3.5.4 shows the coefficients of the method (with the notations of the Butcher tableau). For γ = 1−12√
2
γ γ 0
1−γ 1−2γ γ 1/2 1/2
Table 3.5.4: Coefficients used in the two-stage DIRK method.
the method is L-stable and second-order.
The reference solution for the discretized problem was computed by the Matlab’s ODE45 solver and is plotted in Figure 3.5.2.
In the different operator splittings the difference operator on the right-hand side of the semi-discrete problem was split into the sum D+R, where D contained the discretized diffusion and the inhomogeneous boundary conditions, and R the reaction and source terms. The spatial discretization of the diffusion terms and the big difference in the magnitude of the reaction rates give rise to stiffness, also indicated by the big operator norms kDk=O(103) andkRk=O(104).
Figure 3.5.2: The reference solution of the reaction-diffusion problem.
First we compare the operator splittings of second order accuracy. We emphasize on the fact that here both sub-operators are stiff, while most studies are restricted to the case where one of the operators is stiff, the other is non-stiff. E.g., in [147] it is shown that in that case the ordering of the sub-operators in the sequential and Strang-Marchuk splitting can affect the accuracy considerably, and moreover, always the stiff operator should be put to the end. In our example, however, two stiff operators are present, so both sequences D-R and R-D should be tested for the Richardson-extrapolated sequential splitting and both D-R-D and R-D-R for the Strang-Marchuk splitting. (Obviously, the symmetrically weighted sequential splitting does not require ordering considerations.)
The sub-problems were solved by four different time integration methods: 1) the implicit Euler method, 2) the explicit Euler method, 3) the midpoint method and 4) the two-step DIRK method. During the experiments, the Richardson-extrapolated sequen-tial splitting proved to be unstable for methods 2) and 3). Tables 3.5.5–3.5.6 show the maximum norms of the errors at the end of the time interval for methods 1) and 4).
Note that the time steps were halved in these experiments, therefore a factor around 0.5 (in parentheses) corresponds to first-order convergence, while a factor around 0.25 to second-order convergence.
τ Riseq. D-R Riseq. R-D swss SM
1/10 7.07(-3) 5.86(-4) 4.80(-2) 1.26(-2)
1/20 3.67(-3) (0.52) 1.78(-4) (0.30) 2.40(-2) (0.50) 7.52(-3) (0.60) 1/40 1.86(-3) (0.51) 4.96(-5) (0.28) 1.20e(-2) (0.50) 4.28(-3) (0.57) 1/80 9.21(-4) (0.50) 1.31(-5) (0.26) 6.02(-3) (0.50) 2.35(-3) (0.55) 1/160 4.59(-4) (0.50) 3.22(-6) (0.25) 3.01e(-3) (0.50) 1.26(-3) (0.54) 1/320 2.22(-4) (0.48) 6.52(-7) (0.20) 1.50(-3) (0.50) 6.63(-4) (0.53)
Table 3.5.5: Comparing the errors of the solutions obtained by the Richardson-extrapolated sequential splittings, symmetrically weighted sequential splitting and Strang-Marchuk splitting (D-R-D) in the reaction-diffusion problem (3.5.25)–(3.5.26) for the im-plicit Euler method.
For the implicit Euler method the Richardson-extrapolated sequential splitting shows the expected second-order convergence in the sequence R-D. The errors obtained in the sequence D-R are one magnitude higher, moreover, here we only obtained first-order convergence (order reduction). The worst results were produced by the symmetrically weighted sequential splitting, which behaves as a first-order method, just like the
Strang-Marchuk splitting. (This is not surprising, because a second-order splitting method was combined with a first-order numerical method.) For the Strang-Marchuk splitting we only give the errors for the sequence D-R-D, which generally produced better results.
Table 3.5.6 illustrates that it is not worthwhile combining the Richardson-extrapolated sequential splitting with a second-order numerical method. The theoretically derived consistency order is still only two, furthermore, the obtained errors are bigger than for the first-order implicit Euler method. Note that the order of the symmetrically weighted sequential splitting and Strang-Marchuk splitting did not increase to two. This is the result of order reduction caused by stiffness.
τ Riseq. D-R Riseq. R-D swss SM
1/10 2.13(-2) 3.86(-3) 8.43(-2) 1.91(-2)
1/20 1.06(-2) (0.50) 1.86(-3) (0.48) 3.99(-2) (0.47) 9.48(-3) (0.50) 1/40 4.86(-3) (0.46) 8.95(-4) (0.48) 1.86(-2) (0.47) 4.42(-3) (0.47) 1/80 2.51(-3) (0.52) 3.47(-4) (0.39) 8.55(-3) (0.46) 2.25(-3) (0.51) 1/160 1.10(-3) (0.44) 9.50(-5) (0.27) 3.92(-3) (0.46) 1.01(-3) (0.45) 1/320 3.71(-4) (0.34) 2.04(-5) (0.21) 1.80(-3) (0.46) 3.76(-4) (0.37) 1/640 9.42(-5) (0.25) 8.54(-6) (0.42) 8.45(-4) (0.47) 1.15(-4) (0.31)
Table 3.5.6: Comparing the errors of the solutions obtained by the Richardson-extrapolated sequential splittings, symmetrically weighted sequential splitting and Strang-Marchuk splitting (D-R-D) in the reaction-diffusion problem (3.5.25)–(3.5.26) for the two-step DIRK method.
We have also tested the additive splitting and the iterated splitting on this example.
(The iteration number in the iterative splitting was m= 2.)
The errors were computed in the maximum norm for different values of τ, see Ta-ble 3.5.7. The first column contains the values of the splitting time step. In the second column we can see the errors of the DIRK method without splitting; apparently we have a stable, second-order method. The third column gives the errors of the iterative splitting.
One can see that even if an L-stable numerical method was used, the errors grow ex-tremely fast. This instability can also be observed in Figure 3.5.3. The iterative splitting was also run with preconditioning. As the fourth column shows, this did not stabilize the method. The fifth column contains the errors of the additive splitting, which behaves as a stable, first-order method.
τ DIRK Iterative Precond. iterative Additive
1
10 2.3052(-4) 3.9243 3.5312 0.0843
1
20 5.6682(-5) 5745.9 5191 0.0398
1
40 1.4207(-5) 2.0694(+10) 1.8572(+10) 0.0185 Table 3.5.7: Global errors at T = 1/2 for different values of τ.
Remark 3.5.8 We found that taking into account the intermediate values of the stages in the DIRK method stabilized the results. We applied preconditioning in several different ways, namely: i) by using the additive splitting on the whole interval, ii) by using the result of the additive splitting only at the endpoint of each splitting time interval, iii) by
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
−2
−1 0 1x 1010
x
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
−2
−1 0 1x 106
x
Figure 3.5.3: The solution of the iterative splitting for the two components at T = 1/2 (τ = 1/40).
using the result of the sequential splitting at the endpoint of each splitting time interval.
The best results, shown in Table 3.5.8, were obtained in case ii).
τ Iterative Precond. iterative
1
10 3.2904(-4) 2.9508(-4)
1
20 6.7570(-5) 6.2572(-5)
1
40 1.5476e(-5) 1.4767(-5)
Table 3.5.8: Global errors at T = 1/2 for different values of τ. The preconditioning was done by using the result of the Richardson-extrapolated sequential splitting only at the endpoint of each splitting time interval.