- Email: [email protected]
On the additive splitting procedures and their computer realization
Available online at www.sciencedirect.com
Applied Mathematical Modelling 32 (2008) 1552–1569 www.elsevier.com/locate/apm
On the additive splitting procedures and their computer realization Istva´n Farago´
a,*
, Per Grove Thomsen b, Zahari Zlatev
c
a
b
Eo¨tvo¨s Lora´nd University, Institute of Mathematics, Pa´zma´ny P. s. 1/c, 1117 Budapest, Hungary Institute for Informatics and Mathematical Modelling, Technical University of Denmark, DK-2800 Lyngby, Denmark c National Environmental Research Institute, Frederiksborgvej 399, P.O. Box 358, DK-4000 Roskilde, Denmark Received 1 December 2006; received in revised form 1 April 2007; accepted 12 April 2007 Available online 6 July 2007
Abstract Two additive splitting procedures are defined and studied in this paper. It is shown that these splitting procedures have good stability properties. Some other splitting procedures, which are traditionally used in mathematical models used in many scientific and engineering fields, are sketched. All splitting procedures are tested by using six different numerical methods for solving differential equations. Many conclusions, which are related both to the comparison of the additive splitting procedures with the other splitting procedures and to the influence of the numerical methods for solving differential equations on the accuracy of the splitting procedures, are drawn. 2007 Elsevier Inc. All rights reserved. Keywords: Operator splitting; Numerical methods; Time discretization; Stability
1. Statement of the problem In our investigation we will assume that there are only two operators, i.e. we will demonstrate our methods on the (abstract) Cauchy problem of the form: dwðtÞ ¼ ðA þ BÞwðtÞ; dt
t 2 ð0; T ; wð0Þ ¼ w0 :
ð1Þ
The theoretical results are proved under an assumption that the involved operators are bounded linear operators, and hence the exact solution is wðtÞ ¼ expðtðA þ BÞÞwð0Þ. The experimental results indicate that good results should also be expected when this assumption is relaxed.
*
Corresponding author. Tel./fax: +36 1 209 0555. E-mail address: [email protected] (I. Farago´).
0307-904X/$ - see front matter 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2007.04.017
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1553
2. Traditional operator splitting methods In this section we give a short overview of the different known operator splitting methods which are used in many different applications. For more details, see [1–3]. 2.1. Sequential splitting If the sub-problems involving operators A and B from (1) are treated one after the other, then the resulting algorithm is called a sequential splitting procedure. It is very often worthwhile to describe how the different splitting procedures are to be applied in practice by • giving the order in which the simple operators A and B are applied, and • indicating the splitting time-step size s, which is actually used. In our particular case, see again (1), the application of the sequential splitting procedure at a given splitting time-step can be described by sequence ðAÞs ;
ðBÞs :
ð2Þ
It is necessary to explain how the two sub-models are coupled. Assume that n splitting time-steps have successfully been performed and the next splitting time-step, time-step n + 1, has to be carried out. The approximation obtained at time-step n is used as a starting approximation when the first sub-model is treated. The approximation obtained at the end of computations related to the first sub-model is used as a starting value for the second sub-model. The approximation obtained when the computations related to the second submodel are accomplished is accepted as an approximation of the solution of problem (1) at time-step n + 1. In this way everything is prepared to start the computations related to time-step n + 2. It is necessary to explain how to start the computations at time-step 1, but this is not causing problems, because it is assumed that w(0) = w0 is given; see again (1). It should be noted here that if we change the order of the application of the operators, then the results will normally not be the same, i.e. the sequence (A)s, (B)s is in general different from the sequence (B)s, (A)s. The sequential splitting is in general leading to a numerical approximation of order one. The implication of this fact is that as a rule it is not advisable to use numerical algorithms of order higher than one in the treatment of the sub-problems involving the simpler operators A and B when a sequential splitting procedure is to be used. 2.2. Marchuk–Strang splitting Sometimes it is desirable to apply more accurate splitting procedures. Accuracy of order two can be achieved in the following way. Consider an arbitrary splitting time-step, say step n (i.e. the computations are to be carried out from t = tn to t = tn+1 = tn + s). Assume that the sub-models are treated as follows: • Carry out computations by using the first operator from t = tn to t = tn + 0.5s. • Use the second operator to perform computations from t = tn to t = tn + s. • Perform computations from t = tn + 0.5s to t = tn + s by applying again the first operator. This splitting procedure was proposed in 1968 simultaneously by Marchuk and Strang (see [4,5]). It is also called symmetric splitting. In the notation used in the previous sub-section the calculations at an arbitrary splitting time-step can be described by the sequence: ðAÞ0:5s ;
ðBÞs ;
ðAÞ0:5s :
ð3Þ
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1554
The same idea, as in the previous section, can be applied to couple the used operators: (i) we start with the approximation obtained at time-step n, (ii) the approximation obtained when the computations related to (A)0.5s are accomplished is taken as a starting value for the second sub-model (B)s, (iii) the approximation obtained when the computations related to second (A)0.5s are accomplished is taken as a starting value for the third sub-model and (iv) the approximation obtained when the computations related to (B)s are accomplished is taken as an approximation of the problem (1) at time-step n + 1. It is clear that after (iv) we are ready to proceed by performing the computations at splitting time-step n + 2. As in the previous section it should be noted here that if we change the order of the application of the operators, then the results will normally not be the same, i.e. the sequence (A)0.5s, (B)s(A)0.5s is in general different from the sequence (B)0.5 s, (A)s(B)0.5s. 2.3. Weighted sequential splitting Consider again (1). Assume that the computations at the splitting time-step under consideration are carried out by applying successively the sub-problems as shown by the following two sequences of operators: ðAÞs ;
ðBÞs
ð4Þ
ðBÞs ;
ðAÞs ;
ð5Þ
and
i.e. we perform first a sequential splitting time-step starting with the sub-problem containing operator A and after that proceed with a sequential splitting time-step starting with operator B. Denote by wforward(tn + s) the result obtained when the splitting time-step (4) is finished and by wbackward(tn + s) the corresponding result obtained when the computations involved in (5) are accomplished. Then an approximation: wnþ1 ¼ hwforward ðtn þ sÞ þ ð1 hÞwbackward ðtn þ sÞ
ð6Þ
of the solution w(tn + s) of (1) at t = tn+1 = tn + s is calculated using some weighting parameter h and one can proceed with the calculations for the next splitting time-step. The splitting procedures based on the performance of the calculations at each splitting time-step using the formulae (4)–(6) are called weighted sequential splitting procedures. Very often these splitting procedures are used with h = 0.5. The order of accuracy of the weighted sequential procedures is in general one, but experimental results indicate that these procedures are sometimes more stable than the simple sequential procedure. If h = 0.5, then the order of the weighted sequential splitting procedure is two. The sub-problems involved in (4) and (5) can be coupled in the same way as this was done for the sequential splitting procedure in Section 2. 2.4. Weighted Marchuk–Strang splitting A weighted Marchuk–Strang splitting procedure can be obtained from the ordinary Marchuk–Strang splitting procedure in the same way as the weighted sequential procedure was obtained from the ordinary sequential procedure in the previous sub-section. Assume that the computations at the splitting time-step under consideration are carried out by applying successively the sub-problems as shown by the following two sequences of operators: ðAÞ0:5s ;
ðBÞs ;
ðAÞ0:5s
ð7Þ
ðBÞ0:5s ;
ðAÞs ;
ðBÞ0:5s
ð8Þ
and
i.e. we perform first a Marchuk–Strang splitting time-step starting with the sub-problem containing operator A and after that proceed with a Marchuk–Strang splitting time-step starting with operator B.
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1555
Denote again by wforward(tn + s) the result obtained when the splitting time-step (7) is finished and by wbackward(tn + s) the corresponding result obtained when the computations involved in (8) are accomplished. Then an approximation: wnþ1 ¼ hwforward ðtn þ sÞ þ ð1 hÞwbackward ðtn þ sÞ
ð9Þ
of the solution w(tn + s) of (1) at t = tn+1 = tn + s is calculated using some weighting parameter h and one can proceed with the calculations for the next splitting time-step. The splitting procedures based on the performance of the calculations at each splitting time-step using the formulae (7)–(9) are called weighted Marchuk–Strang splitting procedures or weighted symmetric splitting procedures. Very often these splitting procedures are used with h = 0.5. The order of accuracy of the weighed sequential procedures is in general the same as the order of the ordinary Marchuk–Strang splitting procedures (i.e. two), but experimental results indicate that these procedures are often more stable than the simple sequential procedure. The sub-problems involved in (7) and (8) can be coupled in the same way as this was done for the sequential splitting procedure in Section 2.2. 2.5. Extension for the case of more than two operators The splitting procedures sketched in the previous four sub-sections have been derived for the case where the right-hand-side of the original problems contains as in (1) a sum of two operators. In principle at least, all these procedures can be defined for the case where the right-hand-side of the original problem is a sum of more than two operators (see, for example [1]). 3. Additive splitting This method is based on a simple idea: we solve the different sub-problems by using the same initial function. We obtain the split solution by the use of these results and the initial condition. The algorithm is the following: dwn1 ðtÞ ¼ Awn1 ðtÞ; ðn 1Þs < t 6 ns; dt wn1 ððn 1ÞsÞ ¼ wNsp ððn 1ÞsÞ
ð10Þ
dwn2 ðtÞ ¼ Bwn2 ðtÞ; ðn 1Þs < t 6 ns; dt wn2 ððn 1ÞsÞ ¼ wNsp ððn 1ÞsÞ:
ð11Þ
and
Then the split solution at the mesh-points is defined as wNsp ðnsÞ ¼ wn1 ðnsÞ þ wn2 ðnsÞ wNsp ððn 1ÞsÞ:
ð12Þ
Here n ¼ 1; 2; . . . ; N , where wNsp ð0Þ ¼ w0 . We note that we use a superscript oˆno¨ to indicate the splitting step at which the computations are carried out and a superscript oˆNo¨ to show how many splitting steps are to be performed. One can see that the algorithm can be parallelized on the operator level in a natural way (like the symmetrically weighted sequential splitting) because it uses the same initial element. The local splitting error for bounded linear operators can be investigated directly. (We remark that, by definition, the investigations are performed on the sub-interval [0, s]. However, it is not difficult to see that the same assertions in the following statements remain valid when the sub-interval [0, s] is replaced by any other sub-interval ½ðn 1Þs; ns.)
1556
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
Theorem 1. The additive splitting is a first-order accurate splitting method for bounded linear operators. Proof. The solution of the additive splitting at t = s is defined as wNsp ðsÞ ¼ ½expðAsÞ þ expðBsÞ Iw0 ; where I denotes the identity operator. Hence, we get 1 1 wNsp ðsÞ ¼ I þ As þ A2 s2 þ I þ Bs þ B2 s2 I w0 þ Oðs3 Þ 2 2 1 ¼ I þ ðA þ BÞs þ ðA2 þ B2 Þs2 w0 þ Oðs3 Þ: 2
ð13Þ
ð14Þ
Comparing this expression with the similar Taylor expansion of the exact solution, we get the local splitting error 1 Errsp ðsÞ ¼ ððAB þ BAÞs2 Þw0 þ Oðs3 Þ; 2 which proves the statement.
ð15Þ
h
The following statement shows that the additive splitting approximates the exact solution not only at the mesh-points. Corollary 2. The additive splitting approximates the exact solution on the whole split time-interval [0, s]. Proof. The exact solution of the additive splitting at t 2 ½0; s is wNsp ðtÞ ¼ ½expðAtÞ þ expðBtÞ Iw0 :
ð16Þ
Hence, at any point of the time-interval we have 1 Errsp ðtÞ ¼ ððAB þ BAÞt2 Þw0 þ Oðt3 Þ; 2 which proves the statement.
ð17Þ
h
4. Modified additive splitting method For the additive splitting introduced by (10)–(12) the critical point is the proof of the stability. The problem is based on the definition wNsp ðtÞ in (12): in fact, we should estimate in the norm the operator exp(At) + exp(Bt) I and to show that in norm it is bounded linear by one. Due to the presence of the subtraction in the formula, the triangle inequality is not useful for this aim even in the case when the sub-problems are contractive. Therefore, we modify the method. We execute the separated splitting steps with double operators, and then we compute the splitting approximation as the arithmetical mean of the results. Consequently, the algorithm reads as follows: dwn1 ðtÞ ¼ 2Awn1 ðtÞ; ðn 1Þs < t 6 ns; dt wn1 ððn 1ÞsÞ ¼ wNsp ððn 1ÞsÞ
ð18Þ
dwn2 ðtÞ ¼ 2Bwn2 ðtÞ; ðn 1Þs < t 6 ns; dt wn2 ððn 1ÞsÞ ¼ wNsp ððn 1ÞsÞ:
ð19Þ
and
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1557
Then the split solution at the mesh-points is defined as 1 wNsp ðnsÞ ¼ ðwn1 ðnsÞ þ wn2 ðnsÞÞ: ð20Þ 2 Here, as before n ¼ 1; 2; . . . ; N , while wNsp ð0Þ ¼ w0 . The investigation of the local splitting error (for bounded linear operators) for the modified additive splitting (18)–(20) is similar as it was done in Theorem 1. It should be mentioned that it is possible to replace 2A and 2B with A and B in (18) and (19), respectively. If such a replacement is performed then the integration should be carried out on the interval ðn 1Þs < t 6 ðn þ 1Þs and (20) should be replaced with wNsp ðnsÞ ¼ 0:5ðwn1 ððn þ 1ÞsÞ þ wn2 ððn þ 1ÞsÞÞ. If the operators are linear, then the approach sketched above is equivalent to the approach described by (18)– (20). For non-linear operators the results obtained by using these two approaches will in general be different. Theorem 3. For bounded linear operators the modified additive splitting is a first-order accurate splitting method. Proof. The solution of the modified additive splitting at t = s is defined as 1 wNsp ðsÞ ¼ ½expð2AsÞ þ expð2BsÞw0 : 2 Hence, we get wNsp ðsÞ ¼
ð21Þ
1 1 1 2 2 I þ 2As þ ð2AsÞ þ I þ 2Bs þ ð2BsÞ w0 þ Oðs3 Þ 2 2 2
¼ ðI þ ðA þ BÞs þ ðA2 þ B2 Þs2 Þw0 þ Oðs3 Þ:
ð22Þ
Comparing this expression with the similar Taylor expansion of the exact solution, we get the local splitting error 1 2 Errsp ðsÞ ¼ ððA BÞ s2 Þw0 þ Oðs3 Þ; 2
ð23Þ
which proves the statement. h In the following we analyze the approximation property of the method on the whole time-interval [0, s]. Corollary 4. The additive splitting approximates the exact solution on the whole split time-interval [0, s]. Proof. The exact solution of the modified additive splitting at t 2 [0, s] is 1 wNsp ðtÞ ¼ ½expð2AtÞ þ expð2BtÞw0 : 2 Hence, at any point of the time-interval we have 1 2 Errsp ðtÞ ¼ ðA BÞ t2 w0 þ Oðt3 Þ; 2 which proves the statement. h
ð24Þ
ð25Þ
Remark 5. For the discrete case, a similar approach was suggested in [6]. However, our approach is different because we apply the splitting for the continuous problem. 5. Convergence of the modified additive splitting for contractive bounded linear operators The introduced additive splittings, as any operator splitting method, can be considered as numerical timediscretization methods. Hence, the question of its convergence is a crucial question during the application. Therefore it is necessary to investigate the consistency and stability of the splitting procedure. We assume that the exponentials of the operators tA and tB are contractive, i.e the relations
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1558
k expðtAÞk 6 1 and k expðtBÞk 6 1
ð26Þ
are valid. The notion of consistency is closely related to the local splitting error. Let us denote by Tsp(s) the operator which transforms wNsp ðtn Þ into wNsp ðtnþ1 Þ, i.e. T sp ðsÞwNsp ðtn Þ ¼ wNsp ðtnþ1 Þ:
ð27Þ
Obviously, for the additive splitting T sp ðsÞ ¼ expðsAÞ þ expðsBÞ I;
ð28Þ
while for the modified additive splitting it is 1 T sp ðsÞ ¼ ðexpð2sAÞ þ expð2sBÞÞ: 2 As it is known, a splitting method is called consistent numerical method to the problem (1), when T sp ðsÞ I lim ðA þ BÞ wðtÞ ¼0 s!0 s
ð29Þ
ð30Þ
for all t 2 [0, T] and the convergence is uniform in t. (Here w(t) denotes the solution of the problem (1).) Remark 6. For the solution of the well-posed abstract Cauchy problem (1) the relation wðt þ sÞ wðtÞ lim ðA þ BÞ wðtÞ ¼0 s!0 s
ð31Þ
holds for all t 2 ½0; T and the convergence is uniform in t. Therefore, using (30), a splitting method is consistent if and only if the relation wðt þ sÞ T sp ðsÞwðtÞ ¼0 ð32Þ lim s!0 s for all t 2 ½0; T and the convergence is uniform in t. Substituting t = 0 into (32), we get the relation Errsp ðsÞ ¼0 lim s!0 s
ð33Þ
which coincides with the condition Errsp ðsÞ ¼ Oðspþ1 Þ
ð34Þ
with some p > 0. Consequently, when an operator splitting is consistent then its order can be defined by the rate of the convergence in (30). However, the converse is not true: the relation (34) does not imply automatically the consistency. It is said that a splitting method is stable when the relationship kðT sp ðsÞÞn w0 k < 1
ð35Þ
holds for all ns 2 ½0; T and fixed given initial element w0. For the well-posed problem (1) the Lax theorem states the following: for a consistent numerical method the stability is the necessary and sufficient condition of the convergence. Theorem 7. Under the contractivity condition (26) both the additive and the modified additive splittings are consistent. Proof. First we show the validity of (30) for Tsp(s) defined in (29). We have
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1559
1 ðexpð2sAÞ þ expð2sBÞÞ I T sp ðsÞ I 2 ðA þ BÞ expðtðA þ BÞÞw0 ðA þ BÞ wðtÞ ¼ s s ¼ ½sðA2 þ B2 Þ þ Oðs2 Þ expðtðA þ BÞÞw0 :
ð36Þ
According to (26), kexp(t(A + B))k 6 1. Hence, the right side of (36) converges to zero as s tends to zero, uniformly in t. This results in the consistency of the modified additive splitting.The proof of the statement for the additive splitting is similar. h Due to the consistency property, we have to show the stability. However, as it was already mentioned, for the additive splitting (10)–(12) to do it even for contractive operators, it is a difficult task. But for the modified additive splitting, in the contractive case we can show it directly. Theorem 8. Under the contractivity condition (26) the modified additive splitting is contractive and hence it is stable. Proof. The statement is an obvious consequence of the assumptions and of the formula (20). h Hence, using the Lax theorem, we obtained Theorem 9. Under the contractivity condition (26) the modified additive splitting is convergent. At the end of this section we remark that the definitions of the additive splittings can in an obvious way be extended for non-linear operators, but the results of the theorems will in general not hold when the operators are non-linear. Nevertheless, the numerical results in Section 7 indicate that one should expect reasonably good results also in the case where the operators are non-linear. 6. General discussion of the splitting procedures 6.1. Computing time Assume here that all splitting procedures are used with the same splitting time-stepsize s and that the computing time is the most important factor. Under these assumptions the following conclusions can be drawn. • From a computational point of view, the sequential splitting is the best one (i.e. it is the cheapest one with regard to computing time). • If the computational work for the treatment of the sub-problem involving operator B is much larger than the computational work for treating the sub-problem involving operator A, then the Marchuk–Strang splitting procedure is only slightly more expensive computationally than the sequential splitting procedure. In some environmental models this is precisely the case. In advection-chemistry models (see, for example [7,1]) the chemical part often takes about 95% of the computing time. It is clear that for such models the Marchuk–Strang splitting procedure is quite competitive with the sequential splitting procedure with regard to the computing time used. • The remaining splitting procedures are much more expensive with regards to the computing time than the sequential splitting procedure. The computing time will be increased by a factor at least equal to two when any of these splitting procedures is applied instead of the sequential splitting procedure.
6.2. Accuracy of the results Assume that the accuracy of the results is the most important issue. Then it is perhaps a good idea to apply the Marchuk–Strang splitting procedure or the weighted Marchuk–Strang splitting procedure, because second-order accuracy can be achieved when these splitting procedures are used. It should be emphasized, however, that the use of these two splitting procedures imposes some requirements to the numerical methods that are to be used in the three sub-problems arising when the splitting procedure selected is applied to (1): in
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1560
general the numerical methods for solving differential equations that are to be applied in each sub-problems must be of order at least equal to two. The remaining splitting procedures are of order one. This means, let us re-iterate, that these splitting procedures should not be used when the accuracy required is a critical factor. 6.3. Robustness of the results Some experimental results indicate that the weighted splitting procedures are sometimes more robust than the remaining splitting procedures (i.e. in some difficult situations when the solution changes very quickly these two splitting procedures are giving rather good results; such situations occur during sun-sets and sunrises when some chemical models are run, see [1]). 6.4. Proving theoretical properties It is important to prove some theoretical results (as, for example, convergence results, stability results, etc.). If such results are proved, then it will be clear for which classes of problems it is possible to guarantee that the splitting procedures will ensure good results. Some properties of the additive splitting procedures were proved in this paper. Much more similar results are highly desirable. 6.5. General conclusion The short discussion given in the previous four sub-section shows very clearly that the choice of the best splitting procedure is by far not an easy task. Different requirements arise in the different situations and very often these requirements are working in opposite directions. Therefore, in each particular situation one has first to find out what is the most important requirement and to select the splitting procedure according to this most important requirement. It is also necessary to remember that it is necessary to select a good splitting procedure, but it might turn out that this is not sufficient. The numerical methods for solving differential equations, which have to be used during the numerical treatment of every sub-problem, are also very important and one must be very careful when such methods are to be used in combination with splitting procedures. 7. Numerical results An atmospheric chemical scheme has been used in the numerical experiments. This scheme involves 56 chemical species and it is used in several large air pollution models. For this study, it is important to emphasize that the scheme is described mathematically by a system of non-linear ODEs (more details can be found in [7,1]). The system of ODEs can be written in the following form: dw ¼ f ðt; wÞ; dt
t 2 ½21600; 172800; wð21600Þ ¼ wstart ;
ð37Þ
where wstart is a given vector containing the initial values of the concentrations at the starting point t = 21,600 s. The computations are started at 6:00 in the morning and finishing at 24:00 in the next day (i.e. the length of the time-interval is 42 h). It must be stressed that (a) the system of ODEs (37) is non-autonomous depending indirectly on t through input parameters such as temperature, pressure, humidity, etc.) and (b) the time-interval contains changes from day to night and from night to day (it is important to include such changes because the photo-chemical reactions are deactivated during the night and activated during the day; this causes rapid changes of some components of w). It is assumed that this problem is split into two sub-problems. The first sub-problem contains all chemical reactions in which the ozone compound is participation. The second sub-problem contains all other chemical reactions. The problem was solved without using splitting and by using five splitting procedures: (i) sequential splitting, (ii) symmetric splitting, (iii) weighted sequential splitting, (iv) weighted symmetric splitting and (v) additive splitting. The implementation of the first four splitting procedures is straight-forward. The additive
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1561
splitting procedure has been implemented as follows. Assume that the computations have been completed for some time-step n; ðn ¼ 1; 2; . . .Þ. Then the computations needed to obtain an approximation wn+1 to the exact solution w(tn+1) of (37) are carried out in three steps: (1) Step 1. Carry out two steps with the first sub-problem (proceeding from tn to tn+1) and taking as a starting value in the first of the two steps the accepted approximation wn at t = tn). (2) Step 2. Carry out two steps with the second sub-problem (proceeding from tn to tn+1) and taking again as a starting value in the first of the two steps the accepted approximation w(tn) at t = tn ). (3) Step 3. Calculate the mean value of the two vectors calculated in Step 1 and Step 2 and accept this value as an approximation wn+1 of w(tn+1) at time-step n + 1. Thus, everything is prepared to start the calculations at the next time-step, time-step n + 2. It remains to explain how to start the calculations at the first time-step. This is not causing problems because the initial value of the solution of (37) is given. Six numerical methods for solving systems of ODEs were used in the experiments: (1) (2) (3) (4) (5) (6)
BE, the Backward Euler method (first-order, L-stable; see [8,9]), IMP, the Implicit Mid-point Rule (second-order, A-stable; see [8,9]), RK2, the Two-stage Diagonally Implicit Runge–Kutta Method (second-order, L-stable; see [10]), RK5, the Three-stage Fully Implicit Runge–Kutta Method (fifth-order, L-stable; see [9]), ROS, the Rosenbrock Method (second-order, L-stable; see [11]), TR, the Trapezoidal Rule (second-order, A-stable; see [8,9]),
Numerical results showing the accuracy achieved when the selected splittings procedures are run by using the above six numerical methods for solving ODEs will be presented in this section. The accuracy is measured by ! jwmodel wref s;n s;n j ERROR ¼ max ; ð38Þ 6 n¼0;1;...;168;s¼1;2;...;56 jðjwref s;n j; 10 Þ when the error from all 56 species is measured and ! jwmodel wref s;n s;n j ; ERRORs ¼ max 6 n¼0;1;...;168 jðjwref s;n j; 10 Þ
ð39Þ
when the error for a selected species, species number s, is measured. Some chemical species, as for example the hydroxyl radical, vary in a very wide range. It is shown in Fig. 1 that the concentration of the hydroxyl radical is nearly zero during the night, while it is quickly growing to more than 107 molecules per cubic centimetre during the day. Therefore the test with (38) is much more stringent than the second test (39), when the latter test is applied for species that do not vary too much (such as ozone). 7.1. Accuracy results when no splitting procedure is used We shall start with the case where no splitting is used. Accuracy results obtained by the six methods for solving systems of ODEs, which have been listed above, are given in Tables 1a and 1b. Several important conclusions can be drawn by using the results in Tables 1a and 1b. (1) The Backward Euler Method performs as expected (i.e. increasing the number of time-steps with a factor of two results in an increase of the accuracy with a factor of two). (2) The two A-stable methods (the Implicit Mid-Point Rule and the Trapezoidal Rule) have difficulties when this problem is solved (the expected accuracy is achieved only when the time-stepsize becomes very
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1562
6
OH
x 10
18
16
14
12
10
8
6
4
2
0 5
10
15
20
25
30
35
40
45
50
Fig. 1. The time-variation of the hydroxyl radical.
Table 1a Accuracy obtained when no splitting procedure is used with six numerical methods for solving systems of ODEs NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
3.19E1 1.51E1 7.37E2 3.64E2 1.81E2 9.00E3 4.49E3 2.24E3 1.12E3 5.57E4
3.25E1 5.00E2 8.84E3 1.94E2 2.16E2 2.04E2 1.45E2 3.56E3 1.30E5 5.12E6
6.17E1 2.81E1 1.34E1 6.54E2 3.23E2 1.61E2 8.01E3 4.00E3 2.00E3 9.99E4
6.47E5 3.71E6 3.48E6 3.48E6 3.48E6 3.48E6 3.48E6 3.48E6 3.48E6 3.48E6
1.05E 0 4.37E1 1.99E1 9.48E2 4.62E2 2.29E2 1.14E2 5.67E3 2.83E3 1.41E3
3.10E1 4.78E2 2.38E2 2.37E2 2.32E2 2.12E2 1.49E2 3.67E3 1.34E5 5.14E6
small). This fact indicates that it is desirable to apply L-stable methods for problems arising in atmospheric chemistry (excluding the cases where the time-stepsize is very small). (3) The two second-order L-stable methods (the Diagonally Implicit Runge–Kutta Method and the Rosenbrock-type Method) are in fact performing as first-order methods when the atmospheric chemistry problem is solved. The reason for this behavior is clear when the Rosenbrock-type method is used (this method is originally designed for autonomous problem, while the chemical tests is a non-autonomous ODE system; it is difficult to transform it to a non-autonomous problem but there will still be difficulties because the right-hand function f depends implicitly, through the temperature, on t). It is not very clear why the Diagonally Implicit Runge–Kutta Method performs as a first-order method for this particular problem. (4) The reference solution, which is used to evaluate the error, is of order 106. If the Fifth-order Fully Implicit Runge–Kutta Method is used, then this accuracy is achieved when the time-stepsize is 75 s
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1563
Table 1b The ratios between the current error and the errors obtained in the previous run are given for the case where no splitting procedure is used (the ratios should be approximately equal to 2 for first-order methods, to 4 for second-order methods and so on) NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
– 2.11 2.05 2.02 2.01 2.01 2.00 2.00 2.00 2.01
–
– 2.20 2.10 2.05 2.02 2.01 2.01 2.00 2.00 2.00
– 17.21 1.07 1.00 1.00 1.00 1.00 1.00 1.00 1.00
– 2.40 2.20 2.10 2.05 2.02 2.01 2.01 2.00 2.00
–
6.50 5.66 0.46 0.90 1.06 1.41 4.07 273.85 2.54
6.49 2.01 1.00 1.02 1.09 1.42 4.06 273.88 2.61
(i.e. the number of time-steps is 2048). Therefore, no improvement could be registered when the timestepsize is further reduced. (5) The results indicate that the L-stable Fifth-order Fully Implicit Runge–Kutta Method is probably the best choice when no splitting is used in the solution of problems arising in atmospheric chemistry. 7.2. Accuracy results when the sequential splitting procedure is used Accuracy results obtained when the sequential splitting procedure is used with the same six methods for solving systems of ODEs, as those used in the previous sub-section, are given in Tables 2a and 2b. Several important conclusions can be drawn by using the results in Tables 2a and 2b (1) The Backward Euler Method performs as expected (i.e. increasing the number of time-steps with a factor of two results in an increase of the accuracy with a factor of two). (2) The two A-stable methods (the Implicit Mid-Point Rule and the Trapezoidal Rule) have difficulties when this problem is solved (the accuracy ratios vary in an irregular way). This fact indicates that it is desirable to apply L-stable methods for problems arising in atmospheric chemistry also when a sequential splitting is used. (3) The two second-order L-stable methods (the Diagonally Implicit Runge–Kutta Method and the Rosenbrock-type Method) are in fact performing as first-order methods when the atmospheric chemistry problem is solved. The same reasons as those given in the case where no splitting is used can be given to explain this behavior. (4) The use of the Fifth-order Fully Implicit Runge–Kutta Method does not lead to an increase of the rate of convergence for the combined method (splitting + ODE solver). The combined method is of order one. Table 2a The same as Table 1a, but for the case where the sequential splitting procedure is used NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
3.20E1 1.52E1 7.39E2 3.64E2 1.81E2 9.01E3 4.50E3 2.25E3 1.12E3 5.59E4
3.25E1 5.00E1 1.66E2 1.94E2 2.16E2 2.04E2 1.55E2 3.56E3 2.56E4 1.27E4
6.17E1 2.81E1 1.34E1 6.55E2 3.23E2 1.61E2 8.01E3 4.00E3 2.00E3 9.99E4
8.05E2 3.66E2 1.73E2 8.44E3 4.16E3 2.07E3 1.03E3 5.14E4 2.56E4 1.27E4
1.05E 0 4.37E1 1.99E1 9.48E2 4.63E2 2.29E2 1.14E2 5.67E3 2.83E3 1.41E3
3.10E1 4.78E2 2.38E2 2.37E2 2.32E2 2.12E2 1.49E2 3.36E3 2.56E4 1.27E4
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1564
Table 2b The same as Table 1b, but for the case where the sequential splitting procedure is used NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
– 2.11 2.06 2.03 2.01 2.01 2.00 2.00 2.01 2.00
–
– 2.20 2.10 2.05 2.02 2.01 2.01 2.00 2.00 2.00
– 2.20 2.12 2.05 2.03 2.03 2.01 2.00 2.01 2.02
– 2.40 2.20 2.10 2.05 2.02 2.01 2.01 2.00 2.00
– 6.49 2.01 1.00 1.02 1.09 1.42 4.43 13.13 2.02
6.50 3.01 0.86 0.90 1.06 1.32 4.35 13.91 2.02
(5) The results indicate that the Backward Euler Method is probably the best choice when no splitting is used in the solution of problems arising in atmospheric chemistry. There are no great gains when numerical methods for solving systems of ODEs the order of which is higher than two are used in a combination with the sequential splitting procedure. 7.3. Accuracy results when symmetric splitting procedure is used Accuracy results, which have been obtained by using the symmetric (the Marchuk–Strang) splitting procedure applied with the six methods for solving systems of ODEs used in the previous two sub-sections, are given in Tables 3a and 3b. Also in this case it is possible to draw several important conclusions by using the results that are given in Tables 3a and 3b: (1) The Backward Euler Method, the Diagonally Implicit Runge–Kutta Method and the Rosenbrock-type Method perform as numerical methods of order one (for the first of these methods this should be expected, for the next two methods see the comments after Table 1). (2) The two A-stable methods (the Implicit Mid-Point Rule and the Trapezoidal Rule) have difficulties when this problem is solved (the accuracy ratios vary in an irregular way). This fact indicates that it is more suitable to apply L-stable methods for problems arising in atmospheric chemistry also when a sequential splitting is used. (3) The use of the Fifth-order Fully Implicit Runge–Kutta Method together with the symmetric splitting gives a combined method that performs closely to a method of order two. (4) If the accuracy is the critical parameter, then the results indicate that the Fifth-order L-stable Runge– Kutta Method is in general the best choice (however, the two A-stable methods also give good results when the time-stepsize is sufficiently small). Table 3a The same as Table 1a, but for the case where the symmetric (the Marchuk–Strang) splitting procedure is used NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
1.52E1 7.40E2 3.65E2 1.81E2 9.03E3 4.51E3 2.25E3 1.12E3 5.59E4 2.78E4
5.00E2 8.87E2 1.94E3 2.16E2 2.04E2 1.55E2 3.56E2 1.30E3 5.12E5 3.48E6
2.81E1 1.34E1 6.55E2 3.23E2 1.61E2 8.01E3 4.00E3 2.00E3 9.99E4 4.99E4
3.30E2 1.24E2 5.24E3 1.78E3 5.06E4 1.32E4 3.34E5 8.37E6 3.47E6 3.47E6
4.34E1 1.99E1 9.48E2 4.68E2 2.29E2 1.14E2 5.67E3 2.83E3 1.41E3 7.06E4
4.78E2 2.38E2 2.37E2 2.32E2 2.12E2 1.49E2 3.67E3 1.34E5 5.14E6 3.47E6
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1565
Table 3b The same as Table 1b, but for the case where the symmetric (the Marchuk–Strang) splitting procedure is used NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
– 2.05 2.03 2.02 2.00 2.00 2.00 2.01 2.00 2.01
–
– 2.10 2.05 2.03 2.01 2.01 2.00 2.00 2.00 2.00
– 2.66 2.37 2.94 3.52 3.83 3.95 3.99 2.41 1.00
– 2.20 2.10 2.05 2.02 2.01 2.01 2.00 2.00 2.00
–
0.56 4.57 0.90 1.06 1.32 4.36 274.84 2.54 1.47
2.01 1.00 1.02 1.09 1.42 4.43 273.88 2.61 1.48
(5) It should be pointed out here that in most of the cases the use of the sequential splitting with a stepsize 0.5nt is leading to nearly the same accuracy as that obtained when the symmetric splitting is used with a stepsize nt. It is not very easy to explain this result (probably this is caused by the fact that the major part of the error is due to the errors in the second sub-problem). 7.4. Accuracy results when the weighted sequential splitting procedure is used Accuracy results obtained when the weighted sequential splitting procedure is used with h = 0.5 and with the same six methods for solving systems of ODEs, as those used in the previous sub-section, are given in Tables 4a and 4b. Important conclusions related to the accuracy of the weighted sequential splitting procedure when it is used with different numerical methods for solving systems of ODEs can be drawn by using the results in Tables 4a and 4b: (1) The Backward Euler Method, the Diagonally Implicit Runge–Kutta Method and the Rosenbrock-type Method perform as numerical methods of order one (for the first of these methods this should be expected, for the next two methods see the comments after Table 1). (2) The two A-stable methods (the Implicit Mid-Point Rule and the Trapezoidal Rule) have difficulties when this problem is solved (the accuracy ratios vary in an irregular way). This fact indicates that it is more suitable to apply L-stable methods for problems arising in atmospheric chemistry also when a sequential splitting is used. (3) The use of the Fifth-order Fully Implicit Runge–Kutta Method together with the weighted sequential splitting gives a combined method that performs closely to a method of order two.
Table 4a The same as Table 1a, but for the case where the weighted sequential splitting procedure is used with h = 0.5 NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
3.20E1 1.52E1 7.40E2 3.65E2 1.82E2 9.03E3 4.51E3 2.25E3 1.12E3 5.59E4
3.25E1 5.00E1 8.87E2 1.94E2 2.16E2 2.04E2 1.45E2 3.56E3 1.30E5 5.12E6
6.17E1 2.81E1 1.34E1 6.55E2 3.23E2 1.61E2 8.01E3 4.00E3 2.00E3 9.99E4
6.45E2 2.20E2 7.95E3 3.11E3 1.10E3 3.48E4 1.01E4 2.74E5 7.18E6 3.47E6
1.05E 0 4.37E1 1.99E1 9.48E2 4.68E2 2.29E2 1.14E2 5.67E3 2.81E3 1.41E3
3.10E1 2.78E2 2.38E2 2.37E2 2.32E2 2.12E2 1.49E2 3.51E3 1.34E5 5.14E6
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1566
Table 4b The same as Table 1b, but for the case where the weighted sequential splitting procedure is used with h = 0.5 NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
– 2.11 2.05 2.03 2.02 2.00 2.00 2.00 2.01 2.00
–
– 2.20 2.10 2.05 2.02 2.01 2.01 2.00 2.00 2.00
– 2.93 2.77 2.56 2.83 3.16 3.45 3.69 3.82 2.07
– 2.40 2.20 2.10 2.05 2.02 2.01 2.01 2.00 2.00
–
6.50 3.01 4.52 0.90 1.06 1.32 4.36 273.85 2.02
6.49 2.01 1.00 1.02 1.09 1.42 4.24 261.94 2.61
(4) If the accuracy is the critical parameter, then the results indicate that the Fifth-order L-stable Runge– Kutta Method is in general the best choice (however, the two A-stable methods also give good results when the time-stepsize is sufficiently small). (5) The results obtained when the weighted sequential splitting is used are very similar to those obtained by using the sequential splitting. 7.5. Accuracy results when weighted symmetric splitting procedure is used Accuracy results, which have been obtained by using the weighted symmetric (weighted Marchuk–Strang) splitting procedure applied with h = 0.5 and with the six methods for solving systems of ODEs used in the previous two sub-sections, are given in Tables 5a and 5b. Also in this case it is possible to draw several important conclusions by using the results that are given in Tables 5a and 5b: (1) The Backward Euler Method performs as numerical methods of order one, but the results are much more accurate than the corresponding results in the previous tables (including here the results obtained when no splitting is used). (2) The two A-stable methods (the Implicit Mid-Point Rule and the Trapezoidal Rule) perform as methods of order two (thus, the application of the weighted splitting has some stabilizing effect when A-stable methods are used). (3) The use of the Fifth-order Fully Implicit Runge–Kutta Method together with the weighted sequential splitting gives a combined method that performs closely to a method of order two. (4) The results indicate that the two A-stable second-order methods, the Implicit Mid-point Rule and the Trapezoidal Rule, are the best choice when the weighted symmetric splitting is to be used for this par-
Table 5a The same as Table 1a, but for the case where the symmetric (the Marchuk–Strang) splitting procedure is used with h = 0.5 NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
2.33E1 1.12E2 5.51E2 2.73E2 1.36E2 6.76E3 3.38E3 1.69E3 8.41E4 4.19E4
1.16E1 2.68E2 6.42E3 1.57E3 3.88E4 9.64E5 2.40E5 5.98E6 3.49E6 3.48E6
4.49E1 2.08E1 9.98E2 4.89E2 2.42E2 1.20E2 6.01E3 3.00E3 1.50E3 7.49E4
2.80E2 9.85E3 3.95E3 1.38E3 4.19E4 1.14E4 2.95E5 7.49E6 3.47E6 3.47E6
7.40E1 3.18E1 1.47E1 7.06E2 2.46E2 1.71E2 8.53E3 4.25E3 2.12E3 1.06E3
1.62E2 6.16E3 2.07E3 6.17E4 1.69E4 4.45E5 1.23E5 3.55E6 3.50E6 3.48E6
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1567
Table 5b The same as Table 1b, but for the case where the symmetric (the Marchuk–Strang) splitting procedure is used with h = 0.5 NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
– 2.08 2.03 2.02 2.01 2.01 2.00 2.00 2.01 2.01
– 4.33 4.17 4.09 4.05 4.02 4.02 4.01 1.71 1.00
– 2.16 2.08 2.04 2.02 2.02 2.00 2.00 2.00 2.00
– 2.84 2.49 2.87 3.29 3.67 3.86 3.94 2.16 1.00
– 2.33 2.16 2.08 2.04 2.02 2.00 2.01 2.00 2.00
– 2.63 2.98 3.35 3.65 3.83 3.90 3.18 1.01 1.01
ticular problem. More experiments and/or some theoretical investigations are needed before these two methods can be recommended as a general choice when the weighted symmetric (the weighted Marchuk–Strang) splitting is used with h = 0.5.
7.6. Accuracy results when the additive splitting procedure is used Accuracy results obtained when the additive splitting procedure is used with the same six methods for solving systems of ODEs, as those used in the previous sub-section, are given in Tables 6a and 6b. The particular implementation of the additive splitting procedure, which is used here, was described in the beginning of this section. In a similar way as in the previous sub-sections, important conclusions related to the accuracy of the additive splitting procedure when it is used with different numerical methods for solving systems of ODEs can be drawn by using the results in Tables 6a and 6b: (1) The Backward Euler Method performs as a numerical method of order one, but the results are less accurate than the corresponding results in the previous tables. (2) The two A-stable methods (the Implicit Mid-Point Rule and the Trapezoidal Rule) have great difficulties when this problem is solved (the accuracy ratios vary in a very irregular way). This fact indicates that only L-stable methods should be used for problems arising in atmospheric chemistry when the additive splitting is used. (3) The L-stable methods together with additive splitting give a combined method which is of order one also when the ODE methods are of higher order. This is why it is probably most appropriate to use the additive splitting procedure with first-order numerical methods for solving systems of ODEs although the application of methods of higher order results in slightly more accurate results. Table 6a The same as Table 1a, but for the case where the additive splitting procedure is used NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
1.10E 0 3.25E1 1.53E1 7.43E2 3.66E2 1.82E2 9.04E3 4.51E3 2.25E3 1.12E3
2.35E 1 1.18E 1 5.92E 0 7.70E1 3.78E1 1.73E1 5.36E2 8.34E2 3.08E2 1.32E2
1.88E 0 5.45E1 1.92E1 7.88E2 3.55E2 1.68E2 8.20E3 4.04E3 2.00E3 1.00E3
1.10E 0 2.41E1 8.91E2 4.21E2 2.05E2 1.01E2 5.01E3 2.50E3 1.25E3 6.22E4
2.37E 0 7.10E1 2.58E1 1.08E1 4.95E2 2.36E2 1.16E2 5.71E3 2.84E3 1.42E3
2.15E 1 1.08E 1 5.38E 0 2.68E 0 1.31E 0 1.58E1 1.61E1 8.79E2 1.20E1 1.31E2
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1568
Table 6b The same as Table 1b, but for the case where the additive splitting procedure is used NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
– 3.38 2.12 2.06 2.03 2.01 2.01 2.00 2.00 2.01
– 1.99 1.99 7.69 2.04 2.18 3.23 0.64 2.71 2.33
– 3.45 2.84 2.44 2.22 2.11 2.05 2.03 2.02 2.00
– 4.56 2.70 2.12 2.05 2.03 2.02 2.00 2.00 2.01
– 3.34 2.75 2.39 2.18 2.10 2.03 2.03 2.01 2.00
– 1.99 2.01 2.01 2.05 8.29 0.98 1.83 0.73 9.16
7.7. Accuracy results related to ozone when the additive splitting procedure is used As mentioned in the beginning of this section, the accuracy obtained is much better for species which do not vary too much. This is demonstrated by the numerical results shown in Tables 7a and 7b, where accuracy results related to the ozone concentrations are given. First and foremost it is clearly seen that the results for a smoothly varying chemical species can be much better than the results related to all chemical species (compare the results in Tables 6a and 6b with those in Tables 7a and 7b. Several other conclusions can also be drawn: (1) The combined method behaves as a method of order two for all numerical algorithms for solving systems of ODEs. Table 7a The same as Table 6a but for the case where the accuracy of the ozone concentrations is measured NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
4.83E3 2.43E3 1.21E3 6.04E4 3.02E4 1.51E4 7.54E5 3.72E5 1.89E5 9.42E6
2.43E3 1.21E3 6.05E4 3.02E4 1.51E4 7.55E5 3.78E5 1.89E5 9.43E6 4.71E6
2.54E3 1.24E3 6.16E4 3.06E4 1.52E4 7.58E5 3.78E5 1.89E5 9.44E6 4.71E6
2.43E3 1.21E3 6.05E4 3.02E4 1.51E4 7.55E5 3.78E5 1.89E5 9.43E6 4.71E6
2.59E3 1.26E3 6.20E4 3.07E4 1.53E4 7.59E5 3.79E5 1.89E5 9.44E6 4.71E6
2.44E3 1.21E3 6.05E4 3.02E4 1.51E4 7.55E5 3.78E5 1.89E5 9.43E6 4.71E6
Table 7b The same as Table 6b but for the case where the accuracy of the ozone concentrations is measured NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
– 2.00 1.99 2.00 2.00 2.00 2.00 2.03 1.97 2.01
– 2.01 2.00 2.00 2.00 2.00 2.00 2.00 2.00 2.00
– 2.05 2.01 2.01 2.01 2.01 2.01 2.00 2.00 2.00
– 2.01 2.00 2.00 2.00 2.00 2.00 2.00 2.00 2.00
– 2.06 2.03 2.02 2.01 2.02 2.00 2.01 2.00 2.00
– 2.02 2.00 2.00 2.00 2.00 2.00 2.00 2.00 2.00
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1569
(2) The size of the errors does not depend essentially of the order of numerical algorithm for solving systems of ODEs. This indicates that the error induced by the additive splitting procedure is dominant for this particular compound (it should be noted here that this was not the case for the other splitting procedures (see the previous tables). (3) It is perhaps most appropriate to use the Backward Euler Method with the additive splitting procedure. 7.8. General conclusions • It was assumed that the same method has been used in the solution of the system of differential equations arising in the two sub-problem. This is reasonable when the operators involved in these sub-problems are of the same type. It is clear that it will be much more appropriate to use different methods if the operators are of different types (this will be the case when advection-chemistry problems are treated). • The order of the combined method will in general be equal to min(r, q1, q2) where r is the order of the splitting procedure used, while q1 and q2 are the orders of the numerical methods used to treat the two subproblems. The above formula should be taken into account when the splitting procedure and the numerical methods are selected. • It is important to investigate whether the stabilizing effect of the weighted symmetric (the weighted Marchuk–Strang) splitting procedure can be proved for some sufficiently large class of problems.
References [1] Z. Zlatev, I. Dimov, Computational and numerical challenges in environmental modelling, Studies in Computational Mathematics, vol. 13, Elsevier, 2006. ´ . Havasi, Weighted sequential splittings and their analysis, Comput. Math. Appl. 50 (2005) 1017–1031. [2] P. Csomo´s, I. Farago´, A [3] W. Hundsdorfer, J.G. Verwer, Numerical Solution of Time-dependent Advection–Diffusion-Reaction Equations, Springer, Berlin, 2003. [4] G.I. Marchuk, Some application of splitting-up methods to the solution of mathematical physics problems, Applik. Mat. 12 (2) (1968). [5] G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (1968) 505–517. [6] A.A. Samarsky, P.N. Vabisevich, Additive Schemes for the Problems of Mathematical Physics, Nauka, Moscow, 2002 (in Russian). [7] Z. Zlatev, Computer Treatment of Large Air Pollution Models, Kluwer Academic Publishers., Dordrecht, Boston, London, 1995. [8] J.D. Lambert, Numerical Methods for Ordinary Differential Equations, Wiley, Chichester, New York, Brisbane, Toronto, Singapore, 1991. [9] E. Hairer, G. Wanner, Solving Ordinary Differential Equations, II: Stiff and Differential-algebraic Problems, Springer, Berlin, Heidelberg, New York, London, 1987. [10] Z. Zlatev, Modified diagonally implicit Runge–Kutta methods, SIAM J. Sci. Stat. Comput. 2 (1981) 321–334. [11] J.G. Verwer, I.J. Spee, J.G. Blom, W. Hundsdorfer, A second order Rosenbrock method applied to photochemical dispersion problems, SIAM J. Sci. Comput. 20 (1999) 1456–1480.
Applied Mathematical Modelling 32 (2008) 1552–1569 www.elsevier.com/locate/apm
On the additive splitting procedures and their computer realization Istva´n Farago´
a,*
, Per Grove Thomsen b, Zahari Zlatev
c
a
b
Eo¨tvo¨s Lora´nd University, Institute of Mathematics, Pa´zma´ny P. s. 1/c, 1117 Budapest, Hungary Institute for Informatics and Mathematical Modelling, Technical University of Denmark, DK-2800 Lyngby, Denmark c National Environmental Research Institute, Frederiksborgvej 399, P.O. Box 358, DK-4000 Roskilde, Denmark Received 1 December 2006; received in revised form 1 April 2007; accepted 12 April 2007 Available online 6 July 2007
Abstract Two additive splitting procedures are defined and studied in this paper. It is shown that these splitting procedures have good stability properties. Some other splitting procedures, which are traditionally used in mathematical models used in many scientific and engineering fields, are sketched. All splitting procedures are tested by using six different numerical methods for solving differential equations. Many conclusions, which are related both to the comparison of the additive splitting procedures with the other splitting procedures and to the influence of the numerical methods for solving differential equations on the accuracy of the splitting procedures, are drawn. 2007 Elsevier Inc. All rights reserved. Keywords: Operator splitting; Numerical methods; Time discretization; Stability
1. Statement of the problem In our investigation we will assume that there are only two operators, i.e. we will demonstrate our methods on the (abstract) Cauchy problem of the form: dwðtÞ ¼ ðA þ BÞwðtÞ; dt
t 2 ð0; T ; wð0Þ ¼ w0 :
ð1Þ
The theoretical results are proved under an assumption that the involved operators are bounded linear operators, and hence the exact solution is wðtÞ ¼ expðtðA þ BÞÞwð0Þ. The experimental results indicate that good results should also be expected when this assumption is relaxed.
*
Corresponding author. Tel./fax: +36 1 209 0555. E-mail address: [email protected] (I. Farago´).
0307-904X/$ - see front matter 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2007.04.017
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1553
2. Traditional operator splitting methods In this section we give a short overview of the different known operator splitting methods which are used in many different applications. For more details, see [1–3]. 2.1. Sequential splitting If the sub-problems involving operators A and B from (1) are treated one after the other, then the resulting algorithm is called a sequential splitting procedure. It is very often worthwhile to describe how the different splitting procedures are to be applied in practice by • giving the order in which the simple operators A and B are applied, and • indicating the splitting time-step size s, which is actually used. In our particular case, see again (1), the application of the sequential splitting procedure at a given splitting time-step can be described by sequence ðAÞs ;
ðBÞs :
ð2Þ
It is necessary to explain how the two sub-models are coupled. Assume that n splitting time-steps have successfully been performed and the next splitting time-step, time-step n + 1, has to be carried out. The approximation obtained at time-step n is used as a starting approximation when the first sub-model is treated. The approximation obtained at the end of computations related to the first sub-model is used as a starting value for the second sub-model. The approximation obtained when the computations related to the second submodel are accomplished is accepted as an approximation of the solution of problem (1) at time-step n + 1. In this way everything is prepared to start the computations related to time-step n + 2. It is necessary to explain how to start the computations at time-step 1, but this is not causing problems, because it is assumed that w(0) = w0 is given; see again (1). It should be noted here that if we change the order of the application of the operators, then the results will normally not be the same, i.e. the sequence (A)s, (B)s is in general different from the sequence (B)s, (A)s. The sequential splitting is in general leading to a numerical approximation of order one. The implication of this fact is that as a rule it is not advisable to use numerical algorithms of order higher than one in the treatment of the sub-problems involving the simpler operators A and B when a sequential splitting procedure is to be used. 2.2. Marchuk–Strang splitting Sometimes it is desirable to apply more accurate splitting procedures. Accuracy of order two can be achieved in the following way. Consider an arbitrary splitting time-step, say step n (i.e. the computations are to be carried out from t = tn to t = tn+1 = tn + s). Assume that the sub-models are treated as follows: • Carry out computations by using the first operator from t = tn to t = tn + 0.5s. • Use the second operator to perform computations from t = tn to t = tn + s. • Perform computations from t = tn + 0.5s to t = tn + s by applying again the first operator. This splitting procedure was proposed in 1968 simultaneously by Marchuk and Strang (see [4,5]). It is also called symmetric splitting. In the notation used in the previous sub-section the calculations at an arbitrary splitting time-step can be described by the sequence: ðAÞ0:5s ;
ðBÞs ;
ðAÞ0:5s :
ð3Þ
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1554
The same idea, as in the previous section, can be applied to couple the used operators: (i) we start with the approximation obtained at time-step n, (ii) the approximation obtained when the computations related to (A)0.5s are accomplished is taken as a starting value for the second sub-model (B)s, (iii) the approximation obtained when the computations related to second (A)0.5s are accomplished is taken as a starting value for the third sub-model and (iv) the approximation obtained when the computations related to (B)s are accomplished is taken as an approximation of the problem (1) at time-step n + 1. It is clear that after (iv) we are ready to proceed by performing the computations at splitting time-step n + 2. As in the previous section it should be noted here that if we change the order of the application of the operators, then the results will normally not be the same, i.e. the sequence (A)0.5s, (B)s(A)0.5s is in general different from the sequence (B)0.5 s, (A)s(B)0.5s. 2.3. Weighted sequential splitting Consider again (1). Assume that the computations at the splitting time-step under consideration are carried out by applying successively the sub-problems as shown by the following two sequences of operators: ðAÞs ;
ðBÞs
ð4Þ
ðBÞs ;
ðAÞs ;
ð5Þ
and
i.e. we perform first a sequential splitting time-step starting with the sub-problem containing operator A and after that proceed with a sequential splitting time-step starting with operator B. Denote by wforward(tn + s) the result obtained when the splitting time-step (4) is finished and by wbackward(tn + s) the corresponding result obtained when the computations involved in (5) are accomplished. Then an approximation: wnþ1 ¼ hwforward ðtn þ sÞ þ ð1 hÞwbackward ðtn þ sÞ
ð6Þ
of the solution w(tn + s) of (1) at t = tn+1 = tn + s is calculated using some weighting parameter h and one can proceed with the calculations for the next splitting time-step. The splitting procedures based on the performance of the calculations at each splitting time-step using the formulae (4)–(6) are called weighted sequential splitting procedures. Very often these splitting procedures are used with h = 0.5. The order of accuracy of the weighted sequential procedures is in general one, but experimental results indicate that these procedures are sometimes more stable than the simple sequential procedure. If h = 0.5, then the order of the weighted sequential splitting procedure is two. The sub-problems involved in (4) and (5) can be coupled in the same way as this was done for the sequential splitting procedure in Section 2. 2.4. Weighted Marchuk–Strang splitting A weighted Marchuk–Strang splitting procedure can be obtained from the ordinary Marchuk–Strang splitting procedure in the same way as the weighted sequential procedure was obtained from the ordinary sequential procedure in the previous sub-section. Assume that the computations at the splitting time-step under consideration are carried out by applying successively the sub-problems as shown by the following two sequences of operators: ðAÞ0:5s ;
ðBÞs ;
ðAÞ0:5s
ð7Þ
ðBÞ0:5s ;
ðAÞs ;
ðBÞ0:5s
ð8Þ
and
i.e. we perform first a Marchuk–Strang splitting time-step starting with the sub-problem containing operator A and after that proceed with a Marchuk–Strang splitting time-step starting with operator B.
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1555
Denote again by wforward(tn + s) the result obtained when the splitting time-step (7) is finished and by wbackward(tn + s) the corresponding result obtained when the computations involved in (8) are accomplished. Then an approximation: wnþ1 ¼ hwforward ðtn þ sÞ þ ð1 hÞwbackward ðtn þ sÞ
ð9Þ
of the solution w(tn + s) of (1) at t = tn+1 = tn + s is calculated using some weighting parameter h and one can proceed with the calculations for the next splitting time-step. The splitting procedures based on the performance of the calculations at each splitting time-step using the formulae (7)–(9) are called weighted Marchuk–Strang splitting procedures or weighted symmetric splitting procedures. Very often these splitting procedures are used with h = 0.5. The order of accuracy of the weighed sequential procedures is in general the same as the order of the ordinary Marchuk–Strang splitting procedures (i.e. two), but experimental results indicate that these procedures are often more stable than the simple sequential procedure. The sub-problems involved in (7) and (8) can be coupled in the same way as this was done for the sequential splitting procedure in Section 2.2. 2.5. Extension for the case of more than two operators The splitting procedures sketched in the previous four sub-sections have been derived for the case where the right-hand-side of the original problems contains as in (1) a sum of two operators. In principle at least, all these procedures can be defined for the case where the right-hand-side of the original problem is a sum of more than two operators (see, for example [1]). 3. Additive splitting This method is based on a simple idea: we solve the different sub-problems by using the same initial function. We obtain the split solution by the use of these results and the initial condition. The algorithm is the following: dwn1 ðtÞ ¼ Awn1 ðtÞ; ðn 1Þs < t 6 ns; dt wn1 ððn 1ÞsÞ ¼ wNsp ððn 1ÞsÞ
ð10Þ
dwn2 ðtÞ ¼ Bwn2 ðtÞ; ðn 1Þs < t 6 ns; dt wn2 ððn 1ÞsÞ ¼ wNsp ððn 1ÞsÞ:
ð11Þ
and
Then the split solution at the mesh-points is defined as wNsp ðnsÞ ¼ wn1 ðnsÞ þ wn2 ðnsÞ wNsp ððn 1ÞsÞ:
ð12Þ
Here n ¼ 1; 2; . . . ; N , where wNsp ð0Þ ¼ w0 . We note that we use a superscript oˆno¨ to indicate the splitting step at which the computations are carried out and a superscript oˆNo¨ to show how many splitting steps are to be performed. One can see that the algorithm can be parallelized on the operator level in a natural way (like the symmetrically weighted sequential splitting) because it uses the same initial element. The local splitting error for bounded linear operators can be investigated directly. (We remark that, by definition, the investigations are performed on the sub-interval [0, s]. However, it is not difficult to see that the same assertions in the following statements remain valid when the sub-interval [0, s] is replaced by any other sub-interval ½ðn 1Þs; ns.)
1556
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
Theorem 1. The additive splitting is a first-order accurate splitting method for bounded linear operators. Proof. The solution of the additive splitting at t = s is defined as wNsp ðsÞ ¼ ½expðAsÞ þ expðBsÞ Iw0 ; where I denotes the identity operator. Hence, we get 1 1 wNsp ðsÞ ¼ I þ As þ A2 s2 þ I þ Bs þ B2 s2 I w0 þ Oðs3 Þ 2 2 1 ¼ I þ ðA þ BÞs þ ðA2 þ B2 Þs2 w0 þ Oðs3 Þ: 2
ð13Þ
ð14Þ
Comparing this expression with the similar Taylor expansion of the exact solution, we get the local splitting error 1 Errsp ðsÞ ¼ ððAB þ BAÞs2 Þw0 þ Oðs3 Þ; 2 which proves the statement.
ð15Þ
h
The following statement shows that the additive splitting approximates the exact solution not only at the mesh-points. Corollary 2. The additive splitting approximates the exact solution on the whole split time-interval [0, s]. Proof. The exact solution of the additive splitting at t 2 ½0; s is wNsp ðtÞ ¼ ½expðAtÞ þ expðBtÞ Iw0 :
ð16Þ
Hence, at any point of the time-interval we have 1 Errsp ðtÞ ¼ ððAB þ BAÞt2 Þw0 þ Oðt3 Þ; 2 which proves the statement.
ð17Þ
h
4. Modified additive splitting method For the additive splitting introduced by (10)–(12) the critical point is the proof of the stability. The problem is based on the definition wNsp ðtÞ in (12): in fact, we should estimate in the norm the operator exp(At) + exp(Bt) I and to show that in norm it is bounded linear by one. Due to the presence of the subtraction in the formula, the triangle inequality is not useful for this aim even in the case when the sub-problems are contractive. Therefore, we modify the method. We execute the separated splitting steps with double operators, and then we compute the splitting approximation as the arithmetical mean of the results. Consequently, the algorithm reads as follows: dwn1 ðtÞ ¼ 2Awn1 ðtÞ; ðn 1Þs < t 6 ns; dt wn1 ððn 1ÞsÞ ¼ wNsp ððn 1ÞsÞ
ð18Þ
dwn2 ðtÞ ¼ 2Bwn2 ðtÞ; ðn 1Þs < t 6 ns; dt wn2 ððn 1ÞsÞ ¼ wNsp ððn 1ÞsÞ:
ð19Þ
and
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1557
Then the split solution at the mesh-points is defined as 1 wNsp ðnsÞ ¼ ðwn1 ðnsÞ þ wn2 ðnsÞÞ: ð20Þ 2 Here, as before n ¼ 1; 2; . . . ; N , while wNsp ð0Þ ¼ w0 . The investigation of the local splitting error (for bounded linear operators) for the modified additive splitting (18)–(20) is similar as it was done in Theorem 1. It should be mentioned that it is possible to replace 2A and 2B with A and B in (18) and (19), respectively. If such a replacement is performed then the integration should be carried out on the interval ðn 1Þs < t 6 ðn þ 1Þs and (20) should be replaced with wNsp ðnsÞ ¼ 0:5ðwn1 ððn þ 1ÞsÞ þ wn2 ððn þ 1ÞsÞÞ. If the operators are linear, then the approach sketched above is equivalent to the approach described by (18)– (20). For non-linear operators the results obtained by using these two approaches will in general be different. Theorem 3. For bounded linear operators the modified additive splitting is a first-order accurate splitting method. Proof. The solution of the modified additive splitting at t = s is defined as 1 wNsp ðsÞ ¼ ½expð2AsÞ þ expð2BsÞw0 : 2 Hence, we get wNsp ðsÞ ¼
ð21Þ
1 1 1 2 2 I þ 2As þ ð2AsÞ þ I þ 2Bs þ ð2BsÞ w0 þ Oðs3 Þ 2 2 2
¼ ðI þ ðA þ BÞs þ ðA2 þ B2 Þs2 Þw0 þ Oðs3 Þ:
ð22Þ
Comparing this expression with the similar Taylor expansion of the exact solution, we get the local splitting error 1 2 Errsp ðsÞ ¼ ððA BÞ s2 Þw0 þ Oðs3 Þ; 2
ð23Þ
which proves the statement. h In the following we analyze the approximation property of the method on the whole time-interval [0, s]. Corollary 4. The additive splitting approximates the exact solution on the whole split time-interval [0, s]. Proof. The exact solution of the modified additive splitting at t 2 [0, s] is 1 wNsp ðtÞ ¼ ½expð2AtÞ þ expð2BtÞw0 : 2 Hence, at any point of the time-interval we have 1 2 Errsp ðtÞ ¼ ðA BÞ t2 w0 þ Oðt3 Þ; 2 which proves the statement. h
ð24Þ
ð25Þ
Remark 5. For the discrete case, a similar approach was suggested in [6]. However, our approach is different because we apply the splitting for the continuous problem. 5. Convergence of the modified additive splitting for contractive bounded linear operators The introduced additive splittings, as any operator splitting method, can be considered as numerical timediscretization methods. Hence, the question of its convergence is a crucial question during the application. Therefore it is necessary to investigate the consistency and stability of the splitting procedure. We assume that the exponentials of the operators tA and tB are contractive, i.e the relations
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1558
k expðtAÞk 6 1 and k expðtBÞk 6 1
ð26Þ
are valid. The notion of consistency is closely related to the local splitting error. Let us denote by Tsp(s) the operator which transforms wNsp ðtn Þ into wNsp ðtnþ1 Þ, i.e. T sp ðsÞwNsp ðtn Þ ¼ wNsp ðtnþ1 Þ:
ð27Þ
Obviously, for the additive splitting T sp ðsÞ ¼ expðsAÞ þ expðsBÞ I;
ð28Þ
while for the modified additive splitting it is 1 T sp ðsÞ ¼ ðexpð2sAÞ þ expð2sBÞÞ: 2 As it is known, a splitting method is called consistent numerical method to the problem (1), when T sp ðsÞ I lim ðA þ BÞ wðtÞ ¼0 s!0 s
ð29Þ
ð30Þ
for all t 2 [0, T] and the convergence is uniform in t. (Here w(t) denotes the solution of the problem (1).) Remark 6. For the solution of the well-posed abstract Cauchy problem (1) the relation wðt þ sÞ wðtÞ lim ðA þ BÞ wðtÞ ¼0 s!0 s
ð31Þ
holds for all t 2 ½0; T and the convergence is uniform in t. Therefore, using (30), a splitting method is consistent if and only if the relation wðt þ sÞ T sp ðsÞwðtÞ ¼0 ð32Þ lim s!0 s for all t 2 ½0; T and the convergence is uniform in t. Substituting t = 0 into (32), we get the relation Errsp ðsÞ ¼0 lim s!0 s
ð33Þ
which coincides with the condition Errsp ðsÞ ¼ Oðspþ1 Þ
ð34Þ
with some p > 0. Consequently, when an operator splitting is consistent then its order can be defined by the rate of the convergence in (30). However, the converse is not true: the relation (34) does not imply automatically the consistency. It is said that a splitting method is stable when the relationship kðT sp ðsÞÞn w0 k < 1
ð35Þ
holds for all ns 2 ½0; T and fixed given initial element w0. For the well-posed problem (1) the Lax theorem states the following: for a consistent numerical method the stability is the necessary and sufficient condition of the convergence. Theorem 7. Under the contractivity condition (26) both the additive and the modified additive splittings are consistent. Proof. First we show the validity of (30) for Tsp(s) defined in (29). We have
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1559
1 ðexpð2sAÞ þ expð2sBÞÞ I T sp ðsÞ I 2 ðA þ BÞ expðtðA þ BÞÞw0 ðA þ BÞ wðtÞ ¼ s s ¼ ½sðA2 þ B2 Þ þ Oðs2 Þ expðtðA þ BÞÞw0 :
ð36Þ
According to (26), kexp(t(A + B))k 6 1. Hence, the right side of (36) converges to zero as s tends to zero, uniformly in t. This results in the consistency of the modified additive splitting.The proof of the statement for the additive splitting is similar. h Due to the consistency property, we have to show the stability. However, as it was already mentioned, for the additive splitting (10)–(12) to do it even for contractive operators, it is a difficult task. But for the modified additive splitting, in the contractive case we can show it directly. Theorem 8. Under the contractivity condition (26) the modified additive splitting is contractive and hence it is stable. Proof. The statement is an obvious consequence of the assumptions and of the formula (20). h Hence, using the Lax theorem, we obtained Theorem 9. Under the contractivity condition (26) the modified additive splitting is convergent. At the end of this section we remark that the definitions of the additive splittings can in an obvious way be extended for non-linear operators, but the results of the theorems will in general not hold when the operators are non-linear. Nevertheless, the numerical results in Section 7 indicate that one should expect reasonably good results also in the case where the operators are non-linear. 6. General discussion of the splitting procedures 6.1. Computing time Assume here that all splitting procedures are used with the same splitting time-stepsize s and that the computing time is the most important factor. Under these assumptions the following conclusions can be drawn. • From a computational point of view, the sequential splitting is the best one (i.e. it is the cheapest one with regard to computing time). • If the computational work for the treatment of the sub-problem involving operator B is much larger than the computational work for treating the sub-problem involving operator A, then the Marchuk–Strang splitting procedure is only slightly more expensive computationally than the sequential splitting procedure. In some environmental models this is precisely the case. In advection-chemistry models (see, for example [7,1]) the chemical part often takes about 95% of the computing time. It is clear that for such models the Marchuk–Strang splitting procedure is quite competitive with the sequential splitting procedure with regard to the computing time used. • The remaining splitting procedures are much more expensive with regards to the computing time than the sequential splitting procedure. The computing time will be increased by a factor at least equal to two when any of these splitting procedures is applied instead of the sequential splitting procedure.
6.2. Accuracy of the results Assume that the accuracy of the results is the most important issue. Then it is perhaps a good idea to apply the Marchuk–Strang splitting procedure or the weighted Marchuk–Strang splitting procedure, because second-order accuracy can be achieved when these splitting procedures are used. It should be emphasized, however, that the use of these two splitting procedures imposes some requirements to the numerical methods that are to be used in the three sub-problems arising when the splitting procedure selected is applied to (1): in
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1560
general the numerical methods for solving differential equations that are to be applied in each sub-problems must be of order at least equal to two. The remaining splitting procedures are of order one. This means, let us re-iterate, that these splitting procedures should not be used when the accuracy required is a critical factor. 6.3. Robustness of the results Some experimental results indicate that the weighted splitting procedures are sometimes more robust than the remaining splitting procedures (i.e. in some difficult situations when the solution changes very quickly these two splitting procedures are giving rather good results; such situations occur during sun-sets and sunrises when some chemical models are run, see [1]). 6.4. Proving theoretical properties It is important to prove some theoretical results (as, for example, convergence results, stability results, etc.). If such results are proved, then it will be clear for which classes of problems it is possible to guarantee that the splitting procedures will ensure good results. Some properties of the additive splitting procedures were proved in this paper. Much more similar results are highly desirable. 6.5. General conclusion The short discussion given in the previous four sub-section shows very clearly that the choice of the best splitting procedure is by far not an easy task. Different requirements arise in the different situations and very often these requirements are working in opposite directions. Therefore, in each particular situation one has first to find out what is the most important requirement and to select the splitting procedure according to this most important requirement. It is also necessary to remember that it is necessary to select a good splitting procedure, but it might turn out that this is not sufficient. The numerical methods for solving differential equations, which have to be used during the numerical treatment of every sub-problem, are also very important and one must be very careful when such methods are to be used in combination with splitting procedures. 7. Numerical results An atmospheric chemical scheme has been used in the numerical experiments. This scheme involves 56 chemical species and it is used in several large air pollution models. For this study, it is important to emphasize that the scheme is described mathematically by a system of non-linear ODEs (more details can be found in [7,1]). The system of ODEs can be written in the following form: dw ¼ f ðt; wÞ; dt
t 2 ½21600; 172800; wð21600Þ ¼ wstart ;
ð37Þ
where wstart is a given vector containing the initial values of the concentrations at the starting point t = 21,600 s. The computations are started at 6:00 in the morning and finishing at 24:00 in the next day (i.e. the length of the time-interval is 42 h). It must be stressed that (a) the system of ODEs (37) is non-autonomous depending indirectly on t through input parameters such as temperature, pressure, humidity, etc.) and (b) the time-interval contains changes from day to night and from night to day (it is important to include such changes because the photo-chemical reactions are deactivated during the night and activated during the day; this causes rapid changes of some components of w). It is assumed that this problem is split into two sub-problems. The first sub-problem contains all chemical reactions in which the ozone compound is participation. The second sub-problem contains all other chemical reactions. The problem was solved without using splitting and by using five splitting procedures: (i) sequential splitting, (ii) symmetric splitting, (iii) weighted sequential splitting, (iv) weighted symmetric splitting and (v) additive splitting. The implementation of the first four splitting procedures is straight-forward. The additive
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1561
splitting procedure has been implemented as follows. Assume that the computations have been completed for some time-step n; ðn ¼ 1; 2; . . .Þ. Then the computations needed to obtain an approximation wn+1 to the exact solution w(tn+1) of (37) are carried out in three steps: (1) Step 1. Carry out two steps with the first sub-problem (proceeding from tn to tn+1) and taking as a starting value in the first of the two steps the accepted approximation wn at t = tn). (2) Step 2. Carry out two steps with the second sub-problem (proceeding from tn to tn+1) and taking again as a starting value in the first of the two steps the accepted approximation w(tn) at t = tn ). (3) Step 3. Calculate the mean value of the two vectors calculated in Step 1 and Step 2 and accept this value as an approximation wn+1 of w(tn+1) at time-step n + 1. Thus, everything is prepared to start the calculations at the next time-step, time-step n + 2. It remains to explain how to start the calculations at the first time-step. This is not causing problems because the initial value of the solution of (37) is given. Six numerical methods for solving systems of ODEs were used in the experiments: (1) (2) (3) (4) (5) (6)
BE, the Backward Euler method (first-order, L-stable; see [8,9]), IMP, the Implicit Mid-point Rule (second-order, A-stable; see [8,9]), RK2, the Two-stage Diagonally Implicit Runge–Kutta Method (second-order, L-stable; see [10]), RK5, the Three-stage Fully Implicit Runge–Kutta Method (fifth-order, L-stable; see [9]), ROS, the Rosenbrock Method (second-order, L-stable; see [11]), TR, the Trapezoidal Rule (second-order, A-stable; see [8,9]),
Numerical results showing the accuracy achieved when the selected splittings procedures are run by using the above six numerical methods for solving ODEs will be presented in this section. The accuracy is measured by ! jwmodel wref s;n s;n j ERROR ¼ max ; ð38Þ 6 n¼0;1;...;168;s¼1;2;...;56 jðjwref s;n j; 10 Þ when the error from all 56 species is measured and ! jwmodel wref s;n s;n j ; ERRORs ¼ max 6 n¼0;1;...;168 jðjwref s;n j; 10 Þ
ð39Þ
when the error for a selected species, species number s, is measured. Some chemical species, as for example the hydroxyl radical, vary in a very wide range. It is shown in Fig. 1 that the concentration of the hydroxyl radical is nearly zero during the night, while it is quickly growing to more than 107 molecules per cubic centimetre during the day. Therefore the test with (38) is much more stringent than the second test (39), when the latter test is applied for species that do not vary too much (such as ozone). 7.1. Accuracy results when no splitting procedure is used We shall start with the case where no splitting is used. Accuracy results obtained by the six methods for solving systems of ODEs, which have been listed above, are given in Tables 1a and 1b. Several important conclusions can be drawn by using the results in Tables 1a and 1b. (1) The Backward Euler Method performs as expected (i.e. increasing the number of time-steps with a factor of two results in an increase of the accuracy with a factor of two). (2) The two A-stable methods (the Implicit Mid-Point Rule and the Trapezoidal Rule) have difficulties when this problem is solved (the expected accuracy is achieved only when the time-stepsize becomes very
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1562
6
OH
x 10
18
16
14
12
10
8
6
4
2
0 5
10
15
20
25
30
35
40
45
50
Fig. 1. The time-variation of the hydroxyl radical.
Table 1a Accuracy obtained when no splitting procedure is used with six numerical methods for solving systems of ODEs NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
3.19E1 1.51E1 7.37E2 3.64E2 1.81E2 9.00E3 4.49E3 2.24E3 1.12E3 5.57E4
3.25E1 5.00E2 8.84E3 1.94E2 2.16E2 2.04E2 1.45E2 3.56E3 1.30E5 5.12E6
6.17E1 2.81E1 1.34E1 6.54E2 3.23E2 1.61E2 8.01E3 4.00E3 2.00E3 9.99E4
6.47E5 3.71E6 3.48E6 3.48E6 3.48E6 3.48E6 3.48E6 3.48E6 3.48E6 3.48E6
1.05E 0 4.37E1 1.99E1 9.48E2 4.62E2 2.29E2 1.14E2 5.67E3 2.83E3 1.41E3
3.10E1 4.78E2 2.38E2 2.37E2 2.32E2 2.12E2 1.49E2 3.67E3 1.34E5 5.14E6
small). This fact indicates that it is desirable to apply L-stable methods for problems arising in atmospheric chemistry (excluding the cases where the time-stepsize is very small). (3) The two second-order L-stable methods (the Diagonally Implicit Runge–Kutta Method and the Rosenbrock-type Method) are in fact performing as first-order methods when the atmospheric chemistry problem is solved. The reason for this behavior is clear when the Rosenbrock-type method is used (this method is originally designed for autonomous problem, while the chemical tests is a non-autonomous ODE system; it is difficult to transform it to a non-autonomous problem but there will still be difficulties because the right-hand function f depends implicitly, through the temperature, on t). It is not very clear why the Diagonally Implicit Runge–Kutta Method performs as a first-order method for this particular problem. (4) The reference solution, which is used to evaluate the error, is of order 106. If the Fifth-order Fully Implicit Runge–Kutta Method is used, then this accuracy is achieved when the time-stepsize is 75 s
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1563
Table 1b The ratios between the current error and the errors obtained in the previous run are given for the case where no splitting procedure is used (the ratios should be approximately equal to 2 for first-order methods, to 4 for second-order methods and so on) NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
– 2.11 2.05 2.02 2.01 2.01 2.00 2.00 2.00 2.01
–
– 2.20 2.10 2.05 2.02 2.01 2.01 2.00 2.00 2.00
– 17.21 1.07 1.00 1.00 1.00 1.00 1.00 1.00 1.00
– 2.40 2.20 2.10 2.05 2.02 2.01 2.01 2.00 2.00
–
6.50 5.66 0.46 0.90 1.06 1.41 4.07 273.85 2.54
6.49 2.01 1.00 1.02 1.09 1.42 4.06 273.88 2.61
(i.e. the number of time-steps is 2048). Therefore, no improvement could be registered when the timestepsize is further reduced. (5) The results indicate that the L-stable Fifth-order Fully Implicit Runge–Kutta Method is probably the best choice when no splitting is used in the solution of problems arising in atmospheric chemistry. 7.2. Accuracy results when the sequential splitting procedure is used Accuracy results obtained when the sequential splitting procedure is used with the same six methods for solving systems of ODEs, as those used in the previous sub-section, are given in Tables 2a and 2b. Several important conclusions can be drawn by using the results in Tables 2a and 2b (1) The Backward Euler Method performs as expected (i.e. increasing the number of time-steps with a factor of two results in an increase of the accuracy with a factor of two). (2) The two A-stable methods (the Implicit Mid-Point Rule and the Trapezoidal Rule) have difficulties when this problem is solved (the accuracy ratios vary in an irregular way). This fact indicates that it is desirable to apply L-stable methods for problems arising in atmospheric chemistry also when a sequential splitting is used. (3) The two second-order L-stable methods (the Diagonally Implicit Runge–Kutta Method and the Rosenbrock-type Method) are in fact performing as first-order methods when the atmospheric chemistry problem is solved. The same reasons as those given in the case where no splitting is used can be given to explain this behavior. (4) The use of the Fifth-order Fully Implicit Runge–Kutta Method does not lead to an increase of the rate of convergence for the combined method (splitting + ODE solver). The combined method is of order one. Table 2a The same as Table 1a, but for the case where the sequential splitting procedure is used NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
3.20E1 1.52E1 7.39E2 3.64E2 1.81E2 9.01E3 4.50E3 2.25E3 1.12E3 5.59E4
3.25E1 5.00E1 1.66E2 1.94E2 2.16E2 2.04E2 1.55E2 3.56E3 2.56E4 1.27E4
6.17E1 2.81E1 1.34E1 6.55E2 3.23E2 1.61E2 8.01E3 4.00E3 2.00E3 9.99E4
8.05E2 3.66E2 1.73E2 8.44E3 4.16E3 2.07E3 1.03E3 5.14E4 2.56E4 1.27E4
1.05E 0 4.37E1 1.99E1 9.48E2 4.63E2 2.29E2 1.14E2 5.67E3 2.83E3 1.41E3
3.10E1 4.78E2 2.38E2 2.37E2 2.32E2 2.12E2 1.49E2 3.36E3 2.56E4 1.27E4
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1564
Table 2b The same as Table 1b, but for the case where the sequential splitting procedure is used NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
– 2.11 2.06 2.03 2.01 2.01 2.00 2.00 2.01 2.00
–
– 2.20 2.10 2.05 2.02 2.01 2.01 2.00 2.00 2.00
– 2.20 2.12 2.05 2.03 2.03 2.01 2.00 2.01 2.02
– 2.40 2.20 2.10 2.05 2.02 2.01 2.01 2.00 2.00
– 6.49 2.01 1.00 1.02 1.09 1.42 4.43 13.13 2.02
6.50 3.01 0.86 0.90 1.06 1.32 4.35 13.91 2.02
(5) The results indicate that the Backward Euler Method is probably the best choice when no splitting is used in the solution of problems arising in atmospheric chemistry. There are no great gains when numerical methods for solving systems of ODEs the order of which is higher than two are used in a combination with the sequential splitting procedure. 7.3. Accuracy results when symmetric splitting procedure is used Accuracy results, which have been obtained by using the symmetric (the Marchuk–Strang) splitting procedure applied with the six methods for solving systems of ODEs used in the previous two sub-sections, are given in Tables 3a and 3b. Also in this case it is possible to draw several important conclusions by using the results that are given in Tables 3a and 3b: (1) The Backward Euler Method, the Diagonally Implicit Runge–Kutta Method and the Rosenbrock-type Method perform as numerical methods of order one (for the first of these methods this should be expected, for the next two methods see the comments after Table 1). (2) The two A-stable methods (the Implicit Mid-Point Rule and the Trapezoidal Rule) have difficulties when this problem is solved (the accuracy ratios vary in an irregular way). This fact indicates that it is more suitable to apply L-stable methods for problems arising in atmospheric chemistry also when a sequential splitting is used. (3) The use of the Fifth-order Fully Implicit Runge–Kutta Method together with the symmetric splitting gives a combined method that performs closely to a method of order two. (4) If the accuracy is the critical parameter, then the results indicate that the Fifth-order L-stable Runge– Kutta Method is in general the best choice (however, the two A-stable methods also give good results when the time-stepsize is sufficiently small). Table 3a The same as Table 1a, but for the case where the symmetric (the Marchuk–Strang) splitting procedure is used NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
1.52E1 7.40E2 3.65E2 1.81E2 9.03E3 4.51E3 2.25E3 1.12E3 5.59E4 2.78E4
5.00E2 8.87E2 1.94E3 2.16E2 2.04E2 1.55E2 3.56E2 1.30E3 5.12E5 3.48E6
2.81E1 1.34E1 6.55E2 3.23E2 1.61E2 8.01E3 4.00E3 2.00E3 9.99E4 4.99E4
3.30E2 1.24E2 5.24E3 1.78E3 5.06E4 1.32E4 3.34E5 8.37E6 3.47E6 3.47E6
4.34E1 1.99E1 9.48E2 4.68E2 2.29E2 1.14E2 5.67E3 2.83E3 1.41E3 7.06E4
4.78E2 2.38E2 2.37E2 2.32E2 2.12E2 1.49E2 3.67E3 1.34E5 5.14E6 3.47E6
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1565
Table 3b The same as Table 1b, but for the case where the symmetric (the Marchuk–Strang) splitting procedure is used NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
– 2.05 2.03 2.02 2.00 2.00 2.00 2.01 2.00 2.01
–
– 2.10 2.05 2.03 2.01 2.01 2.00 2.00 2.00 2.00
– 2.66 2.37 2.94 3.52 3.83 3.95 3.99 2.41 1.00
– 2.20 2.10 2.05 2.02 2.01 2.01 2.00 2.00 2.00
–
0.56 4.57 0.90 1.06 1.32 4.36 274.84 2.54 1.47
2.01 1.00 1.02 1.09 1.42 4.43 273.88 2.61 1.48
(5) It should be pointed out here that in most of the cases the use of the sequential splitting with a stepsize 0.5nt is leading to nearly the same accuracy as that obtained when the symmetric splitting is used with a stepsize nt. It is not very easy to explain this result (probably this is caused by the fact that the major part of the error is due to the errors in the second sub-problem). 7.4. Accuracy results when the weighted sequential splitting procedure is used Accuracy results obtained when the weighted sequential splitting procedure is used with h = 0.5 and with the same six methods for solving systems of ODEs, as those used in the previous sub-section, are given in Tables 4a and 4b. Important conclusions related to the accuracy of the weighted sequential splitting procedure when it is used with different numerical methods for solving systems of ODEs can be drawn by using the results in Tables 4a and 4b: (1) The Backward Euler Method, the Diagonally Implicit Runge–Kutta Method and the Rosenbrock-type Method perform as numerical methods of order one (for the first of these methods this should be expected, for the next two methods see the comments after Table 1). (2) The two A-stable methods (the Implicit Mid-Point Rule and the Trapezoidal Rule) have difficulties when this problem is solved (the accuracy ratios vary in an irregular way). This fact indicates that it is more suitable to apply L-stable methods for problems arising in atmospheric chemistry also when a sequential splitting is used. (3) The use of the Fifth-order Fully Implicit Runge–Kutta Method together with the weighted sequential splitting gives a combined method that performs closely to a method of order two.
Table 4a The same as Table 1a, but for the case where the weighted sequential splitting procedure is used with h = 0.5 NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
3.20E1 1.52E1 7.40E2 3.65E2 1.82E2 9.03E3 4.51E3 2.25E3 1.12E3 5.59E4
3.25E1 5.00E1 8.87E2 1.94E2 2.16E2 2.04E2 1.45E2 3.56E3 1.30E5 5.12E6
6.17E1 2.81E1 1.34E1 6.55E2 3.23E2 1.61E2 8.01E3 4.00E3 2.00E3 9.99E4
6.45E2 2.20E2 7.95E3 3.11E3 1.10E3 3.48E4 1.01E4 2.74E5 7.18E6 3.47E6
1.05E 0 4.37E1 1.99E1 9.48E2 4.68E2 2.29E2 1.14E2 5.67E3 2.81E3 1.41E3
3.10E1 2.78E2 2.38E2 2.37E2 2.32E2 2.12E2 1.49E2 3.51E3 1.34E5 5.14E6
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1566
Table 4b The same as Table 1b, but for the case where the weighted sequential splitting procedure is used with h = 0.5 NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
– 2.11 2.05 2.03 2.02 2.00 2.00 2.00 2.01 2.00
–
– 2.20 2.10 2.05 2.02 2.01 2.01 2.00 2.00 2.00
– 2.93 2.77 2.56 2.83 3.16 3.45 3.69 3.82 2.07
– 2.40 2.20 2.10 2.05 2.02 2.01 2.01 2.00 2.00
–
6.50 3.01 4.52 0.90 1.06 1.32 4.36 273.85 2.02
6.49 2.01 1.00 1.02 1.09 1.42 4.24 261.94 2.61
(4) If the accuracy is the critical parameter, then the results indicate that the Fifth-order L-stable Runge– Kutta Method is in general the best choice (however, the two A-stable methods also give good results when the time-stepsize is sufficiently small). (5) The results obtained when the weighted sequential splitting is used are very similar to those obtained by using the sequential splitting. 7.5. Accuracy results when weighted symmetric splitting procedure is used Accuracy results, which have been obtained by using the weighted symmetric (weighted Marchuk–Strang) splitting procedure applied with h = 0.5 and with the six methods for solving systems of ODEs used in the previous two sub-sections, are given in Tables 5a and 5b. Also in this case it is possible to draw several important conclusions by using the results that are given in Tables 5a and 5b: (1) The Backward Euler Method performs as numerical methods of order one, but the results are much more accurate than the corresponding results in the previous tables (including here the results obtained when no splitting is used). (2) The two A-stable methods (the Implicit Mid-Point Rule and the Trapezoidal Rule) perform as methods of order two (thus, the application of the weighted splitting has some stabilizing effect when A-stable methods are used). (3) The use of the Fifth-order Fully Implicit Runge–Kutta Method together with the weighted sequential splitting gives a combined method that performs closely to a method of order two. (4) The results indicate that the two A-stable second-order methods, the Implicit Mid-point Rule and the Trapezoidal Rule, are the best choice when the weighted symmetric splitting is to be used for this par-
Table 5a The same as Table 1a, but for the case where the symmetric (the Marchuk–Strang) splitting procedure is used with h = 0.5 NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
2.33E1 1.12E2 5.51E2 2.73E2 1.36E2 6.76E3 3.38E3 1.69E3 8.41E4 4.19E4
1.16E1 2.68E2 6.42E3 1.57E3 3.88E4 9.64E5 2.40E5 5.98E6 3.49E6 3.48E6
4.49E1 2.08E1 9.98E2 4.89E2 2.42E2 1.20E2 6.01E3 3.00E3 1.50E3 7.49E4
2.80E2 9.85E3 3.95E3 1.38E3 4.19E4 1.14E4 2.95E5 7.49E6 3.47E6 3.47E6
7.40E1 3.18E1 1.47E1 7.06E2 2.46E2 1.71E2 8.53E3 4.25E3 2.12E3 1.06E3
1.62E2 6.16E3 2.07E3 6.17E4 1.69E4 4.45E5 1.23E5 3.55E6 3.50E6 3.48E6
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1567
Table 5b The same as Table 1b, but for the case where the symmetric (the Marchuk–Strang) splitting procedure is used with h = 0.5 NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
– 2.08 2.03 2.02 2.01 2.01 2.00 2.00 2.01 2.01
– 4.33 4.17 4.09 4.05 4.02 4.02 4.01 1.71 1.00
– 2.16 2.08 2.04 2.02 2.02 2.00 2.00 2.00 2.00
– 2.84 2.49 2.87 3.29 3.67 3.86 3.94 2.16 1.00
– 2.33 2.16 2.08 2.04 2.02 2.00 2.01 2.00 2.00
– 2.63 2.98 3.35 3.65 3.83 3.90 3.18 1.01 1.01
ticular problem. More experiments and/or some theoretical investigations are needed before these two methods can be recommended as a general choice when the weighted symmetric (the weighted Marchuk–Strang) splitting is used with h = 0.5.
7.6. Accuracy results when the additive splitting procedure is used Accuracy results obtained when the additive splitting procedure is used with the same six methods for solving systems of ODEs, as those used in the previous sub-section, are given in Tables 6a and 6b. The particular implementation of the additive splitting procedure, which is used here, was described in the beginning of this section. In a similar way as in the previous sub-sections, important conclusions related to the accuracy of the additive splitting procedure when it is used with different numerical methods for solving systems of ODEs can be drawn by using the results in Tables 6a and 6b: (1) The Backward Euler Method performs as a numerical method of order one, but the results are less accurate than the corresponding results in the previous tables. (2) The two A-stable methods (the Implicit Mid-Point Rule and the Trapezoidal Rule) have great difficulties when this problem is solved (the accuracy ratios vary in a very irregular way). This fact indicates that only L-stable methods should be used for problems arising in atmospheric chemistry when the additive splitting is used. (3) The L-stable methods together with additive splitting give a combined method which is of order one also when the ODE methods are of higher order. This is why it is probably most appropriate to use the additive splitting procedure with first-order numerical methods for solving systems of ODEs although the application of methods of higher order results in slightly more accurate results. Table 6a The same as Table 1a, but for the case where the additive splitting procedure is used NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
1.10E 0 3.25E1 1.53E1 7.43E2 3.66E2 1.82E2 9.04E3 4.51E3 2.25E3 1.12E3
2.35E 1 1.18E 1 5.92E 0 7.70E1 3.78E1 1.73E1 5.36E2 8.34E2 3.08E2 1.32E2
1.88E 0 5.45E1 1.92E1 7.88E2 3.55E2 1.68E2 8.20E3 4.04E3 2.00E3 1.00E3
1.10E 0 2.41E1 8.91E2 4.21E2 2.05E2 1.01E2 5.01E3 2.50E3 1.25E3 6.22E4
2.37E 0 7.10E1 2.58E1 1.08E1 4.95E2 2.36E2 1.16E2 5.71E3 2.84E3 1.42E3
2.15E 1 1.08E 1 5.38E 0 2.68E 0 1.31E 0 1.58E1 1.61E1 8.79E2 1.20E1 1.31E2
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1568
Table 6b The same as Table 1b, but for the case where the additive splitting procedure is used NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
– 3.38 2.12 2.06 2.03 2.01 2.01 2.00 2.00 2.01
– 1.99 1.99 7.69 2.04 2.18 3.23 0.64 2.71 2.33
– 3.45 2.84 2.44 2.22 2.11 2.05 2.03 2.02 2.00
– 4.56 2.70 2.12 2.05 2.03 2.02 2.00 2.00 2.01
– 3.34 2.75 2.39 2.18 2.10 2.03 2.03 2.01 2.00
– 1.99 2.01 2.01 2.05 8.29 0.98 1.83 0.73 9.16
7.7. Accuracy results related to ozone when the additive splitting procedure is used As mentioned in the beginning of this section, the accuracy obtained is much better for species which do not vary too much. This is demonstrated by the numerical results shown in Tables 7a and 7b, where accuracy results related to the ozone concentrations are given. First and foremost it is clearly seen that the results for a smoothly varying chemical species can be much better than the results related to all chemical species (compare the results in Tables 6a and 6b with those in Tables 7a and 7b. Several other conclusions can also be drawn: (1) The combined method behaves as a method of order two for all numerical algorithms for solving systems of ODEs. Table 7a The same as Table 6a but for the case where the accuracy of the ozone concentrations is measured NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
4.83E3 2.43E3 1.21E3 6.04E4 3.02E4 1.51E4 7.54E5 3.72E5 1.89E5 9.42E6
2.43E3 1.21E3 6.05E4 3.02E4 1.51E4 7.55E5 3.78E5 1.89E5 9.43E6 4.71E6
2.54E3 1.24E3 6.16E4 3.06E4 1.52E4 7.58E5 3.78E5 1.89E5 9.44E6 4.71E6
2.43E3 1.21E3 6.05E4 3.02E4 1.51E4 7.55E5 3.78E5 1.89E5 9.43E6 4.71E6
2.59E3 1.26E3 6.20E4 3.07E4 1.53E4 7.59E5 3.79E5 1.89E5 9.44E6 4.71E6
2.44E3 1.21E3 6.05E4 3.02E4 1.51E4 7.55E5 3.78E5 1.89E5 9.43E6 4.71E6
Table 7b The same as Table 6b but for the case where the accuracy of the ozone concentrations is measured NT
BE
IMP
RK2
RK5
ROS
TR
1008 2016 4032 8064 16,128 32,256 64,512 129,024 258,048 512,096
– 2.00 1.99 2.00 2.00 2.00 2.00 2.03 1.97 2.01
– 2.01 2.00 2.00 2.00 2.00 2.00 2.00 2.00 2.00
– 2.05 2.01 2.01 2.01 2.01 2.01 2.00 2.00 2.00
– 2.01 2.00 2.00 2.00 2.00 2.00 2.00 2.00 2.00
– 2.06 2.03 2.02 2.01 2.02 2.00 2.01 2.00 2.00
– 2.02 2.00 2.00 2.00 2.00 2.00 2.00 2.00 2.00
I. Farago´ et al. / Applied Mathematical Modelling 32 (2008) 1552–1569
1569
(2) The size of the errors does not depend essentially of the order of numerical algorithm for solving systems of ODEs. This indicates that the error induced by the additive splitting procedure is dominant for this particular compound (it should be noted here that this was not the case for the other splitting procedures (see the previous tables). (3) It is perhaps most appropriate to use the Backward Euler Method with the additive splitting procedure. 7.8. General conclusions • It was assumed that the same method has been used in the solution of the system of differential equations arising in the two sub-problem. This is reasonable when the operators involved in these sub-problems are of the same type. It is clear that it will be much more appropriate to use different methods if the operators are of different types (this will be the case when advection-chemistry problems are treated). • The order of the combined method will in general be equal to min(r, q1, q2) where r is the order of the splitting procedure used, while q1 and q2 are the orders of the numerical methods used to treat the two subproblems. The above formula should be taken into account when the splitting procedure and the numerical methods are selected. • It is important to investigate whether the stabilizing effect of the weighted symmetric (the weighted Marchuk–Strang) splitting procedure can be proved for some sufficiently large class of problems.
References [1] Z. Zlatev, I. Dimov, Computational and numerical challenges in environmental modelling, Studies in Computational Mathematics, vol. 13, Elsevier, 2006. ´ . Havasi, Weighted sequential splittings and their analysis, Comput. Math. Appl. 50 (2005) 1017–1031. [2] P. Csomo´s, I. Farago´, A [3] W. Hundsdorfer, J.G. Verwer, Numerical Solution of Time-dependent Advection–Diffusion-Reaction Equations, Springer, Berlin, 2003. [4] G.I. Marchuk, Some application of splitting-up methods to the solution of mathematical physics problems, Applik. Mat. 12 (2) (1968). [5] G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (1968) 505–517. [6] A.A. Samarsky, P.N. Vabisevich, Additive Schemes for the Problems of Mathematical Physics, Nauka, Moscow, 2002 (in Russian). [7] Z. Zlatev, Computer Treatment of Large Air Pollution Models, Kluwer Academic Publishers., Dordrecht, Boston, London, 1995. [8] J.D. Lambert, Numerical Methods for Ordinary Differential Equations, Wiley, Chichester, New York, Brisbane, Toronto, Singapore, 1991. [9] E. Hairer, G. Wanner, Solving Ordinary Differential Equations, II: Stiff and Differential-algebraic Problems, Springer, Berlin, Heidelberg, New York, London, 1987. [10] Z. Zlatev, Modified diagonally implicit Runge–Kutta methods, SIAM J. Sci. Stat. Comput. 2 (1981) 321–334. [11] J.G. Verwer, I.J. Spee, J.G. Blom, W. Hundsdorfer, A second order Rosenbrock method applied to photochemical dispersion problems, SIAM J. Sci. Comput. 20 (1999) 1456–1480.