- Email: [email protected]
Alternating Order Algorithm Based on Stages of the Ceschino Method
February 18 - 20, 2015. Vienna University of Technology, Vienna, 8th Vienna International Conference on Mathematical Modelling Austria 8th Vienna Conference on Mathematical Available onlineModelling at www.sciencedirect.com 8th Vienna18International International Conference on Mathematical Modelling February 20, 2015. Vienna University of Technology, February 18 18 -- 20, 20, 2015. 2015. Vienna Vienna University University of of Technology, Technology, Vienna, Vienna, February Vienna, Austria Austria Austria
ScienceDirect
IFAC-PapersOnLine 48-1 (2015) 198–203 Alternating Order Algorithm Based on Stages of the Ceschino Method Alternating Order Algorithm Based on Stages of the Ceschino Method Alternating Based Stages of Alternating Order Order Algorithm Algorithm Based on on Stages of the the Ceschino Ceschino Method Method Eugeny A. Novikov*, Anton E. Novikov**
Eugeny A. Novikov*, Anton E. Novikov** Eugeny Eugeny A. A. Novikov*, Novikov*, Anton Anton E. E. Novikov** Novikov** * Department of Computational Mathematics, Institute of Computational Modeling of SB RAS, Krasnoyarsk, Russia ** Department of Computational Mathematics, Department of Computational Mathematics, (e-mail: [email protected]). * Department of Computational Mathematics, Institute of of Computational Modeling Modeling of of SB RAS, RAS, Krasnoyarsk, Krasnoyarsk, Russia Russia Institute **Department of Mathematical SupportofforSB Devices andRussia Systems, Institute of Computational Computational SBDiscrete RAS, Krasnoyarsk, (e-mail: Modeling [email protected]). (e-mail: [email protected]). Siberian Federal University, Krasnoyarsk, Russia (e-mail: [email protected]). **Department of Mathematical Support for Discrete Devices and Systems, **Department of Support for Discrete and (e-mail: [email protected]) **Department of Mathematical Mathematical Support for Discrete Devices Devices and Systems, Systems, Siberian Federal University, Krasnoyarsk, Russia Siberian University, Krasnoyarsk, Russia Siberian Federal Federal University, Krasnoyarsk, Russia (e-mail: [email protected]) (e-mail: (e-mail: [email protected]) [email protected]) Abstract: An inequality for stability control for the Ceschino method of the second order is constructed. The First order method with an extended stability domain based on of stages of the order Ceschino scheme is Abstract: An inequality for stability control for the Ceschino method the second is constructed. Abstract: An inequality for stability control for the Ceschino method of the second order is constructed. constructed. Calculations results of stiff problems modeling, confirming the increased efficiency of Abstract: An inequality for stability control for the Ceschino method of the second order is constructed. The First order method with an extended stability domain based on stages of the Ceschino scheme is The First order method with an extended stability domain based on stages of the Ceschino scheme is numerical calculations due to using alternating order are given. The First order method with an extended stability domain based on stages of the Ceschino scheme is constructed. Calculations results of stiff problems modeling, confirming the increased efficiency of constructed. Calculations Calculations results results of of stiff stiff problems problems modeling, modeling, confirming confirming the the increased increased efficiency efficiency of of constructed. numerical calculations due to to using alternating alternating ordercontrol, are given. given. © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: explicit methods, and stability stiff problems. numerical due using order numerical calculations calculations due toaccuracy using alternating order are are given. Keywords: explicit methods, accuracy and stability control, stiff problems. Keywords: Keywords: explicit explicit methods, methods, accuracy accuracy and and stability stability control, control, stiff stiff problems. problems. Jacobi matrix and, consequently, to the increase of 1. INTRODUCTION computational second approach is based on the Jacobi matrix matrix costs. and, The consequently, to the the increase of Jacobi and, consequently, to increase of 1. INTRODUCTION Jacobi matrix and, consequently, to the increase of λ of the Jacobi estimation of the maximum eigenvalue The necessity of applying algorithms based on explicit 1. INTRODUCTION max computational costs. The second approach is based on the 1. INTRODUCTION computational costs. second approach is the computational costs. The The secondusing approach is based based on on methods appears inapplying numerical large-scale stiffonproblems matrix by the power iterations the λincrements of the the the Jacobi estimation of the the maximum eigenvalue The necessity of algorithms based explicit max of λ of the Jacobi estimation of maximum eigenvalue The necessity of applying algorithms based on explicit max λmax of the Jacobi estimation of the maximum eigenvalue solving as shown in [1–3]. Integration methods on a base of The necessity of applying algorithms based on explicit right-hand part of a system of ordinary differential equations methods appears in numerical large-scale stiff problems matrix by the power iterations using the increments of the methods in numerical large-scale stiff problems matrix the iterations using the increments of the implicit orappears semi-implicit numerical schemes usually use the methods appears in[1–3]. numerical large-scale stiff problems matrix bysubsequent the power power iterations using the increments of [7]. the with theby control of ordinary the inequality h ⋅ λmax equations ≤D solving as shown in Integration methods on a base of right-hand part of a system of differential solving as shown in [1–3]. Integration methods on a base of right-hand part of a system of ordinary differential equations decomposition of the Jacobi matrix [2], that is a separate solving as shown in [1–3]. Integration methods on a base of right-hand part of a system of ordinary differential equations estimation doescontrol not leadoftothe theinequality increase ofh computational implicit or semi-implicit numerical schemes usually use the This with the subsequent ⋅λ ≤ D implicit or schemes usually use with subsequent control of h D [7]. [7]. time-consuming task. In numerical suchmatrix a situation, the algorithms implicit or semi-implicit semi-implicit numerical schemes usually use the the costs max ≤ with the the subsequent control of the the inequality inequality h ⋅⋅ λλmax in all the considered situations [3–5, 7]. max ≤ D [7]. decomposition of the Jacobi [2], that is a separate decomposition of the Jacobi matrix [2], that is a separate This estimation does not lead to the increase of computational based on explicit formulas are more efficient if stiffness of decomposition oftask. the Jacobi matrix [2], that is algorithms a separate This estimation does to of time-consuming In such aa situation, This estimation does not not lead lead to the the increase increase of computational computational time-consuming task. In an such situation,tothe the algorithms costs in all the considered situations [3–5, 7]. the problem allows to get approximation a solution in a time-consuming task. In such a situation, the algorithms In this paper, the stability control inequality costs in all the considered situations [3–5, 7]. based on on explicit explicit formulas formulas are are more more efficient efficient if if stiffness stiffness of of costs in all the considered situations [3–5, 7]. for the second based reasonable time [3]. based on explicit formulas are more efficient stiffnessinofa order Ceschino method [8] is constructed using the approach the problem allows to get an approximation to aaifsolution In the stability control inequality for the the problem allows to get an approximation to solution in In this thisinpaper, paper, thefirst stability control inequality forstages the second second the problem allows to get an approximation to a solution in aa In given [7]. The order method based on the of the this paper, the stability inequality for second reasonable time [3]. order Ceschino method [8] control is constructed using thetheapproach approach The step control algorithm is usually based on accuracy reasonable time [3]. order Ceschino method [8] is constructed using the reasonable time [3]. numerical Ceschino formula with a stability interval extended order Ceschino method [8] is constructed using the approach in [7]. The first order method based on the stages of the control of control a numerical scheme. It is natural, because the given given [7]. The first order based the of to 32 in along the real axis constructed. Ainterval newstages integration The step algorithm is usually based on accuracy given in [7]. The firstformula orderismethod method based on on the stages of the the The step control algorithm is usually based on accuracy numerical Ceschino with a stability extended accuracy the approximation the main of numerical The stepofof control algorithm is is usually basedrequirement on accuracy Ceschino formula with aa stability interval extended algorithm of alternating order and step is formulated. The control a numerical scheme. It is natural, because the numerical Ceschino formula with stability interval extended control of numerical scheme. It is natural, because the to 32 along the real axis is constructed. A new integration calculations. However, applying algorithms control ofofaa the numerical scheme.isthe Itthe isintegration natural, because the to 32 the real is constructed. A new integration numerical results of axis stifforder problems solving confirming the accuracy approximation main requirement of to 32 along alongof the real axis is constructed. A formulated. new integration accuracy of the approximation is the main requirement of algorithm alternating and step is The based on explicit methods to stiff problems leads to the loss accuracy of the approximation is the main requirement of algorithm of alternating order and step is formulated. The increase of the efficiency due to alternating order and the calculations. However, applying the integration algorithms algorithm of alternating order and step is formulated. The calculations. However, applying the integration algorithms results of stiff problems solving confirming the of efficiency and reliability [4–5]. The leads reason is loss the numerical calculations. However, applying the integration algorithms numerical results of stiff problems solving confirming the stability control are given. based on explicit methods to stiff problems to the numerical results of stiff problems solving confirming the based on on explicit explicit methods to stiff stiff problems problems leadsonto to the loss increase of the efficiency due to alternating order and the contradiction between the accuracy and stability settling based methods to of efficiency due of efficiency efficiency and reliability [4–5]. The leads reasona the is loss the increase increase of the the are efficiency due to to alternating alternating order order and and the the of and reliability [4–5]. The reason is the stability control given. region which leads to a large amount of the recomputed of efficiency and reliability [4–5]. The reason is the stability control control are are given. contradiction between the accuracy and stability on aa settling stability 2. given. CESCHINO METHOD contradiction between the accuracy and stability on settling solutions, whereas thetothe step is much lessstability thanthethe maximum contradiction between and on a settling region which leads aa accuracy large amount of recomputed region which leads to large amount of the recomputed 2. CESCHINO CESCHINO METHOD allowable one. The additional stability control of a numerical region which leads to a large amount of the recomputed We consider the Cauchy problemMETHOD 2. solutions, whereas the step is much less than the maximum 2. CESCHINO METHOD solutions, whereas the step is much less than the maximum scheme allows to avoid this problem. At the present moment, solutions, whereas the step is much less than the maximum allowable one. The additional stability control of a numerical We consider the Cauchy problem y ′ = f (t, y ) , the y (t0 ) = y0 , problem t0 ≤ t ≤ t k , (1) allowable one. The additional stability of aa numerical there are two common approaches topresent stability We allowable one. The additional stabilityAtcontrol control ofthe numerical We consider consider the Cauchy Cauchy problem scheme allows to avoid this problem. the moment, scheme allows to avoid this problem. At the present moment, control [5–7], obtained applying the Dahlquist test equation. scheme allows to avoid this problem. At thetopresent moment, yy ′′ = f (ty, yand ) , y (ft ) are = yy0 N , t-dimensional ≤ tt ≤ t , (1) where real vector functions, there are two common approaches the stability (1) there are two common approaches to the stability y′ = = ff ((tt,, yy )) ,, yy ((tt000 )) = = y00 ,, tt000 ≤ ≤t ≤ ≤ ttkkk ,, (1) there are two common approaches to the stability control [5–7], obtained applying the Dahlquist test equation. t is an independent variable, which varies over the interval where and y f are N -dimensional real vector functions, The first approach refers with the estimating the maximum control [5–7], [5–7], obtained obtained applying applying the the Dahlquist Dahlquist test test equation. equation. where and y f are N -dimensional real vector functions, control and f are N -dimensional real vector functions, eigenvalue of the refers Jacobi with matrix its norm with the where For the solution of (1) we use theover explicit Runge[ttt0isis, tkan ] . yindependent variable, which varies the interval The first first approach approach the using estimating the maximum maximum an independent The refers with the estimating the t is an independent variable, variable, which which varies varies over over the the interval interval subsequent (in addition to the accuracy control) control of the The first approach refers with the estimating the maximum Kutta eigenvalue of the Jacobi matrix using norm with the . For the the solution of of (1) we we use the the explicit RungeRunge[[ttt00 ,,, tttkk ]]formulas eigenvalue of the matrix using its its with the eigenvalue of faddition the Jacobi Jacobi matrix its norm norm with the For the solution solution of (1) (1) we use use the explicit explicit Runge[ 0 k ] .. For ancontrol) integration stepof and inequality h isusing h(in subsequent to the accuracy control the y ≤ D , where subsequent (in addition to the accuracy control) control of the Kutta formulas yn +1 =formulas yn + pm1k1 + pm 2 k2 + pm 3 k3 + pm 4 k4 , subsequent (in addition to the accuracy control) control of the Kutta Kutta formulas ,, where is an integration step and inequality h h ff yof≤ D D is a norm ODE system Jacobi matrix [6]. Here, is f where is an integration step and inequality h h ≤ D y y inequality h f y ≤ D , where h is an integration step and y = yyn + pm1kk1 + pm 2 kk2⎛+ p 1k + pm 41kk4 ,, ⎞ 2 k2 + m 4 4k kyy1nnn +++=111 = hf y(nntn+ t p+mm 33 kh33 ,+y p+ , (2) = +, yp pnmm)11k, 11 k+ +2 p p=mmhf 2 2⎜+n pm 3 k3 + npm 4 k 4 ,1 ⎟ a ffpositive constant correlated with the size[6]. of Here, a stability D is aa norm of ODE system Jacobi matrix is y D is norm of ODE system Jacobi matrix [6]. Here, is 4 4 ⎞⎠ ⎝ 1 1 y ⎛ D is a norm of ODE system Jacobi matrix [6]. Here, is f 1 1 y ⎛ ⎞ domain. It is = hf hf ( t , y ) , k =hf hf ⎛⎜ ttn + + 1 h,, yyn + + 1 kk1 ⎞⎟ ,, (2) kkk11 = (2) positive constant constant correlated correlated with with the the size size of of aa stability stability = hf ((⎛ttnn ,, yynn1)) ,, kk22 = =hf h, ynn + 4 1 ⎜ t⎞n + 4 h (2) aaawell-known, positive 4 k11 ⎟⎟⎠⎠ , that explicit numerical schemes involve positive constant correlated with the size do of not a stability k13 =hf ⎜ tnn +n h, y2n + k⎜⎝⎝⎝2 ⎟n , 4 4 4 domain. It is ⎠ domain. It is 2 2 ⎠⎞ the Jacobi that matrix. Therefore, thenotapproach domain. It applying is ⎝⎛ 1 well-known, explicit numerical schemes do involve 1 1 ⎛⎜ tn + 1 kk3 = hf ,, yyn + well-known, that explicit numerical schemes do not involve 1h 1 kk2 ⎞⎞⎟ ,, ⎛ = hf t + h + mentioned above leads to the additional calculation of the well-known, that explicit numerical schemes do not involve k33 =hf ⎜⎝⎜ tnn + 2 h, ynn + 2 k22 ⎟⎠⎟ , the Jacobi matrix. Therefore, applying the approach 2 2 the matrix. Therefore, applying the ⎝⎝ 2 2 ⎠⎠ the Jacobi Jacobiabove matrix. Therefore, applying the approach approach mentioned leads to the additional calculation of the mentioned above leads mentioned© 2015, aboveIFAC leads to to the the additional additional calculation calculation of of the the 198 Copyright 2405-8963 © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Peer review©under of International Federation of Automatic Control. Copyright 2015, responsibility IFAC 198 Copyright © 198 10.1016/j.ifacol.2015.05.080 Copyright © 2015, 2015, IFAC IFAC 198
MATHMOD 2015 February 18 - 20, 2015. Vienna, Austria Eugeny A. Novikov et al. / IFAC-PapersOnLine 48-1 (2015) 198–203
k= hf ( tn + h, yn + k1 − 2k2 + 2k3 ) , 4
where X = h ⋅ A . It is easy to see, that the following relations
where h is an integration step, ki , 1 ≤ i ≤ 4 , are the stages of the method, pmi , 1 ≤ i ≤ 4 , are the numerical coefficients, m is the accuracy order of the method. On the coefficients p21 = 1 , p22 = −2 , p23 = 2 , p24 = 0 ,
199
k1 − 2k2 + k3=
(3)
the scheme (2), (3) has the second accuracy order [8].
1 3 X ⋅ yn , 8
1 1 2 k1 ) X ⋅ yn , ( k2 −= 2 8
hold. The estimation of the maximum eigenvalue of the Jacobi matrix of the system (1) can be calculated by the power iterations [3]. We introduce the notation
3. ACCURACY CONTROL The scheme (2) has the fourth order on the coefficients 1 2 1 , p 42 =0 , p 43 = , p 44 = . The estimation of the local 6 3 6 error δn ,2 of the form p 41 =
⎪⎧ ( k1 − 2k2 + k3 )i ⎪⎫ vn = 2 ⋅ max ⎨ ⎬. 1≤ i ≤ N ⎪⎩ ( k2 − k1 )i ⎪⎭
δn ,2 = ( p41 − p21 ) k1 + ( p42 − p22 ) k2 +
For stability control of the Ceschino method the inequality vn ≤ D can be applied, where the value of D bounds the stability interval.
+ ( p43 − p23 ) k3 + ( p44 − p24 ) k4 can be used for the accuracy control. As a result, for the accuracy control we use the inequality
(4)
The stability properties of the Runge-Kutta methods are usually investigated on the scalar test equation y ′ = λy , where λ is some arbitrary complex number, Re ( λ ) < 0 . λ
δn ,2 ≤ ε ,
pretends to be an eigenvalue of the Jacobi matrix of the problem (1). Applying (2), (3) to the solution of y ′ = λy we
where ⋅ is some norm in R N , ε is the required accuracy of calculations. Here, the relation δn ,2 = O ( h 3 ) holds.
get that the stability function Q2 ( x ) of the second order method has the form
In this paper the term "accuracy step" is referred to the step chosen according to the accuracy requirements and the term "stability step" denotes the step, which is chosen according to the stability requirements. We choose the accuracy step hac using the formula h ac= q ⋅ h , where q is determined from
Q2 ( x ) =1 + x +
1 2 1 3 x + x , x= h ⋅ λ , 2 4
and the stability function Q4 ( x ) of the fourth order method
the equation q 3 ⋅ δn ,2 = ε . If q < 1 , we recompute the
has the following form
solution with the step h equal to q ⋅ h . Otherwise, we compute and accept the solution and the predicted step hn +1 is calculated using the formula hn +1= q ⋅ h . The inequality
Q4 ( x ) =1 + x +
1 2 1 3 1 4 x + x + x , x= h ⋅ λ . 2 6 24
The length of the stability interval of the second order method is equal to 2 (fig. 1), whereas for the fourth order scheme it is equal to 2.8 (fig. 2). Therefore, in the inequality vn ≤ D we let D = 2 . Taking into account, that vn = O ( h ) ,
δn ,2 ≤ ε is well-behaved in solving a lot of practical
problems, therefore, we use it in this paper. In this paper the algorithm of alternating step based on the scheme (2), (3) with the inequality for accuracy control δn ,2 ≤ ε will be
the stability step h st can be determined from the formula hst = r ⋅ h , where r is calculated from the equation , D 2. r= ⋅ vn D=
referred to as CESCH42 .
There are three reasons why estimation (4) is rather rough: 1) the maximum eigenvalue may be very close to other eigenvalues; 2) there are only a few iterations in the power iterations method; 3) as the estimation is obtained on the linear problem y ′ = Ay , non-linearity of the problem (1) leads to the additional error in the estimation of the maximum eigenvalue. Since the estimation is rough the stability control is used to limit stepsize. As a result, the predicted stepsize is calculated using the formula
4. STABILITY CONTROL of NUMERICAL SCHEME Now we construct the inequality for stability control of the numerical scheme (2). For this purpose, we apply (2) for the solution of the linear problem y ′ = Ay with some constant matrix A . The first three stages k1 , k2 and k3 applied to this problem have the following form
1 2⎞ ⎛ k= X ⋅ yn , k 2 = 1 ⎜ X + X ⎟ ⋅ yn , 4 ⎝ ⎠ 1 1 ⎛ ⎞ k3 =⎜ X + X 2 + X 3 ⎟ ⋅ yn , 2 8 ⎝ ⎠
{
}
hn +1 = max hn , min ⎡⎣ h ac , h st ⎤⎦ ,
199
(5)
MATHMOD 2015 February 18 - 20, 2015. Vienna, Austria Eugeny A. Novikov et al. / IFAC-PapersOnLine 48-1 (2015) 198–203 200
settling region is much higher than the required one. It is natural, because on this region the previous errors are suppressed due to the stability control, whereas new errors are small due to low values of the derivatives of the solution. In this situation, it is more efficient to perform calculations using some low order method with wider stability domain.
where hn is the latest successful integration step. Notice, that formula (5) is applied to predict the stepsize of hn +1 after successful calculation of a solution with the previous step hn and, therefore, does not lead to the increase of computational costs. If the stability step is less than the latest successful one, it will not be decreased, because the demand of decreasing the stepsize may be caused by rather roughness of the estimation. However, the stepsize will also not be increased, because the possibility of instability of a numerical scheme is not excluded. If the stability step has to be decreased, the latest successful step hn is applied. As a result, the formula (5) is used to choose the stepsize. The given formula allows to stabilize the step behavior on the settling region of the solution, where the stability plays a defining role. It is the existence of the settling area, that significantly restricts the capabilities of the explicit methods on solving stiff problems.
5. FIRST ORDER METHOD Now we construct the first order method yn +1,1 = yn + p11k1 + p12 k 2 + p13 k3 + p14 k 4 ,
(6)
based on the stages of the numerical scheme (2) with a wider stability domain. For this purpose, we apply (6) to the solution of the test equation y ′ = λy . We get
yn +1,1 = Q1 ( x ) yn , where x = h ⋅ λ and the stability function Q1 ( x ) has the form Q1 ( x ) = 1 + ( p11 + p12 + p13 + p14 ) x + 1 ⎛1 ⎞ + ⎜ p12 + p13 + p14 ⎟ x 2 + 2 ⎝4 ⎠ 1 1 ⎛1 ⎞ + ⎜ p13 + p14 ⎟ x3 + p14 x 4 . 8 2 4 ⎝ ⎠
(7)
The requirement of the first accuracy order of the numerical formula (7) leads to the relation p11 + p12 + p13 + p14 = 1.
We use the other coefficients p1i , 1 ≤ i ≤ 4 , to make the stability domain wider. The stability condition for the scheme (6) has the form Q1 ( x ) ≤ 1 . In order to construct a method
Fig. 1. The stability domain of the second order method.
with the maximum stability interval we consider the Chebyshev polynomial
T4 ( z ) = 8 z 4 − 8 z 2 + 1 . It is known, that the polynom T4 ( z ) is less deflected from zero for z ∈ [ −1,1] . We introduce the notation z= 1−
2 x, γ
so that the interval [ γ, 0] is mapped to the interval [ −1,1] . As a result, the polynom T4 ( x ) has the form T4 ( x ) = 1−
32 160 256 128 x + 2 x 2 − 3 x3 + 4 x 4 . γ γ γ γ
It is easy to show [3], that for the polynom T4 ( x )
Fig. 2. The stability domain of the fourth order method.
(8) the
inequality T4 ( x ) ≤ 1 is satisfied in the maximum interval
In this paper, the algorithm of alternating step with the additional stability control of the numerical scheme will be referred to as CESCH42st . The given algorithm is based on a numerical formula of low (the second) accuracy order, therefore, it is aimed at the solution of non-stiff and moderate stiff problems with low accuracy of calculations (of 1% order and lower). From the calculations with the algorithm CESCH42st it follows, that the practical accuracy on a
[ γ, 0 ] ,
γ = −32 (fig. 3). Comparing the relations (7) and (8)
on γ = −32 , we get the coefficients p1i , 1 ≤ i ≤ 4 , of the first order method (2) with the maximum stability interval, i.e. p11 = 200
895 257 31 1 , p12 = , p13 = , p14 = . 2048 512 512 2048
(9)
MATHMOD 2015 February 18 - 20, 2015. Vienna, Austria Eugeny A. Novikov et al. / IFAC-PapersOnLine 48-1 (2015) 198–203
In the fig. 3 the stability domain consists of four parts, which intersect in extreme points of the Chebyshev polynomial. For those problems, which Jacobi matrix eigenvalues have nonzero imaginary part, it may lead to decrease of the calculation efficiency. This effect may be overcome by applying stability polynomials which are close to optimal polynomials. In [3] there is given an algorithm for construction stability polynomials with stability domain extended along the implicit axis.
201
order method are accompanied by (in addition to accuracy control) the control of the inequality vn ≤ 32 , whereas the stepsize is chosen using the formula (5). In this paper the algorithm of alternating order and step is referred to as CESCH42vp . 8. DIFFERENTIAL EQUATIONS of CHEMICAL KINETICS The kinetic scheme of any chemical reaction involves the following elementary reactions (stages) α1,1 x1 + ... + α NR,1 xNR → β1,1 x1 + ... + βNR,1 xNR ,
……… α1,NS x1 + ... + α NR,NS xNR → β1,NS x1 + ... + βNR,NS xNR ,
(10)
where xi , 1 ≤ i ≤ NR are chemical reagents, NR and NS are the number of the reagents and the stages in the reaction, respectively; αij and βij , 1 ≤ i ≤ NR , 1 ≤ j ≤ NS are the stoichiometric coefficients. For each elementary reaction the corresponding stage velocity constant k j , 1 ≤ j ≤ NS is
Fig. 3. The stability domain of the first order method.
given. For the process (10) under the lumped model of the constant-volume isothermal reactor the corresponding ordinary differential equations system has the following form
The stability domain of the first order method (6), (9) along the real axis is 16 times wider than the stability domain of the numerical scheme (2), (3). The computational costs of the first and the second order method are the same. Therefore, for the problems, in which the stepsize is limited mostly due to the stability issues the assumed theoretic increase of the efficiency is equal to 16 times.
C ′ = AT V , C ( 0 ) = C0 . Here, AT is a stoichiometric matrix, C and V are vectors of reagents concentrations and stage velocities, respectively. In case of the reaction flowing in non-isothermal conditions this system also involves the heat balance equation
6. ACCURACY AND STABILITY CONTROL of FIRST ORDER METHOD
T′ =
In the accuracy control inequality we use the following estimation of the local truncation error
QT V − α (T − T01 ) CVT C
,
where T is the mixture temperature in the reactor, T01 is the
δn ,1 = ( p41 − p11 ) k1 + ( p42 − p12 ) k2 +
temperature of the reactor walls, QT is the vector of the
+ ( p43 − p13 ) k3 + ( p44 − p14 ) k4 .
calorific capacities of the stages, CVT is a vector of the heat % / r , α% is the heat-conduction capacities of the reagents, α = αs coefficient, s and r are square of area and reactor volume, the vectors respectively. The superscript T for QT and CVT denotes the transposition. The heat capacities of the reagents and the heat-conduction coefficient may be functions of the reagents concentrations ci , 1 ≤ i ≤ NR , α may depend on the temperature.
For the accuracy control of the numerical formula (6), (9) we can apply the inequality δn ,1 ≤ ε . Since the stability interval of the numerical scheme (6), (9) is limited by 32, we apply for stability control the inequality vn ≤ 32 where vn is determined using the formula (4). 7. ALGORITHM of ALTERNATING ORDER
If the reaction proceeds in a constant-volume isothermal reactor with the substance exchange (an open system, a perfect-mixing reactor) a system of ordinary differential equations has the form
The first order method with the extended stability domain is efficient on the settling regions, where stepsize is limited due to the stability issues. On the transition region it is more efficient to apply the method (2), (3). We can increase the efficiency by applying each method on the region, where it is the most efficient. As a criteria for switching between the methods we can use the inequality for stability control. On calculations with the method (2), (3) we switch to the numerical scheme (6), (9) on the violation of the inequality vn ≤ 2 . On calculations with the first order method we switch back when vn ≤ 2 satisfies. The calculations by the first
C ′ =AT V +
1 ( C p − C ) , C ( 0 ) = C0 . Θ
where C p is a vector of incoming reagent concentrations, Θ is the stay period of the mixture in the reactor, Θ = r / u , u is the mixture space velocity through a reactor. If the
201
MATHMOD 2015 February 18 - 20, 2015. Vienna, Austria Eugeny A. Novikov et al. / IFAC-PapersOnLine 48-1 (2015) 198–203 202
10. NUMERICAL MODELING OF ETHANE PYROLYSIS
reaction proceeds in non-isothermal conditions this system also involves the heat balance equation
T′
QT V − α (T − T01 ) CVT C
The ethane pyrolysis in absence of oxygen is determined by a small sequence of elementary reactions. The mechanism of the ethane pyrolysis has been discussed in scientific literature many times. Here, we use the reaction scheme, which was proposed and investigated in [10]
1 − (T − T02 ) , Θ
where T02 is the incoming mixture temperature in the reactor. The temperature of the reaction mixture can be described by a function of time t and concentrations of reagents ci ,
C2 H 6 → CH 3 +CH3 , CH 3 +C2 H 6 → CH 4 +C 2 H 5 , C2 H 5 → C2 H 4 +H , H+C 2 H 6 → H 2 +C2 H 5 ,
1 ≤ i ≤ NR , i.e. T = T ( t , C ) .
C 2 H 5 +C2 H 5 → C 4 H10 .
9. ALGORITHM FOR GENERATING CHEMICAL KINETIC EQUATIONS
Here, the stage velocity constants have the form = k1 1.34 ⋅ 10−5 , k 2 = 3.73 ⋅10 2 ,= k3 3.69 ⋅103 ,
If an elementary reaction is balancing, its velocity Ws is
k 4 = 3.66 ⋅ 105 и k5 = 1.62 ⋅107 .
equal to a difference of its straight Ws+ and backward Ws−
We denote reagent concentrations in the following way
reactions, i.e. W = Ws+ − Ws− , 1 ≤ s ≤ NS . If there is some s third particle in a reaction, the velocity Vs is calculated using the formulas [9]
c1 = [ C2 H 6 ] , c2 = [ CH3 ] , c3 = [ CH 4 ] , c4 = [ C2 H5 ] ,
Vs = PW s s , Ps =
NR+NI
∑ i =1
c5 = [ C2 H 4 ] , c6 = [ H ] , c7 = [ H 2 ] , c8 = [ C4 H10 ] . Then, the corresponding system involves 8 ordinary differential equations of the form C ′ = AT V , i.e.
εsi ci , 1 ≤ s ≤ NS ,
c1′ = −k1c1 − k2 c1c2 − k4 c1c6 = , c2′ 2k1c1 − k2 c1c2 ,
where εsi , 1 ≤ s ≤ NS , NR + 1 ≤ i ≤ NR+NI , are the efficiencies of the third particles, NI is the number of inert substances, εsi and ci are the efficiencies and concentrations of the inert substances, respectively. The values of components of the vector Ws are determined from the scheme of chemical reaction (10) using the relations + s
W = ks
NR+NI
Πc i =1
αij
− s
, W = k− s
NR+NI
Πc i =1
βij
c3′ = k2 c1c2 , c4′ = k 2 c1c2 − k3 c4 + k 4 c1c6 − 2k5 c42 ,
(11)
c5′ = k3 c4 , c6′ = k3c4 − k4 c1c6 , c7′ = k4 c1c6 , c8′ = k c . 2 5 4
The initial concentration of the ethane c1 = [ C2 H 6 ] is equal to 0.14, other reagent concentrations are equal to zero. The calculations were performed with the accuracy ε = 10−2 . The efficiency of the integration algorithms was estimated by the total number if of calculations of the right-hand part of the problem (11) over an integration interval. The numerical approximation was performed over the interval [0,0.26]. The given problem satisfies the 'classical' definition of stiffness. In the beginning of the interval a transition region (centiseconds) is observed and after there is a settling region (fig. 4).
,
where k s and k− s , 1 ≤ s ≤ NS , are straight and backward velocity stage constants, respectively. Velocity constants of the stages are calculated using the formulas ⎛ Ej ⎞ n = k j Aj T j exp ⎜ − ⎟, ⎝ RT ⎠ where T is the temperature of the mixture in a reactor; Aj , n j and E j / R are the given constants. Notice that, in
general, the velocity constants are not constant in case of a non-isothermal reactor – they depend on temperatures. However, historically, an isothermal reactor was considered earlier than the non-isothermal one and thus k j , 1 ≤ j ≤ NS , at the present moment, are still called constants. The stoichiometric matrix AT with the elements aij is formed from the kinetic scheme (10) conforming to the following rule. The number of a stage is aligned with the column number and the reagent number is aligned with the row number of the matrix AT . If xi serves as an initial reagent, then aij = − αij , if xi is a product, then aij = βij . If xi is both an initial reagent and a product, then aij = −αij + βij . Usually a few amount of reagents are reacting in a stage, i.e. a stoichiometric matrix is sparse.
Fig. 4. Time-dependence of CH 3 concentration.
202
MATHMOD 2015 February 18 - 20, 2015. Vienna, Austria Eugeny A. Novikov et al. / IFAC-PapersOnLine 48-1 (2015) 198–203
The comparism of the efficiencies was performed between the constructed algorithms and the known Merson method (MERSON, [11]). In case of applying explicit methods for stiff problems solving the integration stepsize is limited due to the stability requirements almost over the whole integration interval. As a result, the practical accuracy of a solution is much higher than the defined one. In fig. 4 the graphs of a solution are coincide for the Merson method and for the constructed algorithm CESCH42.
203
6. Shampine, L.M. (1982). Implementation of Rosenbrock methods. ACM Transaction on Mathematical Software, No. 5, 93–113. 7. Novikov, V.A. and Novikov, E.A. (1984) Control of the stability of explicit one - step methods of integrating ordinary differential equations. Soviet Math. Dokl., vol. 30, No. 1, 211–215. 8. Ceschino, F., and Kuntzman, G. (1966) Numerical solution of initial value problems. Prentice-Hall, Englewood Clis, New Jersey. 9. Novikov, E.A. (2010) Numerical modeling of a modified oregonator by the (2,1)-method for solving stiff problems. Numerical methods and programming, vol. 11, No. 2, 123–130 (in Russian). 10. Kulich, D.M., and Taylor, J.E. (1975) Mathematical simulation of the oxygen ethane reaction. Chem. Kinetic, vol. 8, 89–97. 11. Merson, R.H. (1957) An operational methods for integration processes. Proc. of Symp. on Data Processing, Australia, 329–331. 12. Enright, W.H., and Hull, T.E. (1975) Comparing numerical methods for the solutions of systems of ODE’s. BIT, No. 15, 10–48.
For CESCH42 without the stability control we have if = 22 853 , for the algorithm CESCH42st with stability control we have if =20 403 , for the algorithm of alternating order and step CESCH42vp the amount of calculations of the right-hand part is if = 2 588 and for the Merson method it is equal to if = 25 796 . 11. CONCLUSIONS From the numerical results we may draw the following important conclusions. Firstly, the constructed integration algorithm of the second order with accuracy control of calculations and the stability control of the numerical scheme as well as the algorithm of alternating order and step can be applied for the solution of rather stiff problems. Secondly, by the calculation costs the algorithm of alternating order and step CESCH42vp is more efficient than the Merson method almost 10 times, that is due to the stability control of the numerical scheme and calculations with alternating order. It seems that at rather high dimension of the problem (1) the method CESCH42vp may compete with implicit methods at solving moderate stiff problems, because it does not require to decompose the Jacobi matrix. The algorithm CESCH42vp is more efficient at the solution of the twelve test problems [2] and ten examples [12]. ACKNOWLEDGMENT This work was partially supported by Project 14-11-00147 of Russian Scientific Foundation. REFERENCES 1. Hairer, E., Norsett, S.P., and Wanner, G. (1987). Solving ordinary differential equations. Stiff and differentialalgebraic problems, Springer-Verlag, Berlin. 2. Hairer, E., and Wanner, G. (1996). Solving ordinary differential equations. Non stiff problems, SpringerVerlag, Berlin. 3. Novikov, E.A. (1997). Explicit methods for stiff systems, Nauka, Novosibirsk (in Russian). 4. Novikov, E.A., and Shornikov, Yu.V. (2013). Computer simulation of stiff hybrid systems,: NSTU Publishers, Novosibirsk (in Russian). 5. Novikov, E.A., and Zacharov, A.A. (2011). Harmonization of the stability regions in the three-stage explicit RungeKutta. Vestnik TyumGU, No. 7, 187–192 (in Russian).
203
ScienceDirect
IFAC-PapersOnLine 48-1 (2015) 198–203 Alternating Order Algorithm Based on Stages of the Ceschino Method Alternating Order Algorithm Based on Stages of the Ceschino Method Alternating Based Stages of Alternating Order Order Algorithm Algorithm Based on on Stages of the the Ceschino Ceschino Method Method Eugeny A. Novikov*, Anton E. Novikov**
Eugeny A. Novikov*, Anton E. Novikov** Eugeny Eugeny A. A. Novikov*, Novikov*, Anton Anton E. E. Novikov** Novikov** * Department of Computational Mathematics, Institute of Computational Modeling of SB RAS, Krasnoyarsk, Russia ** Department of Computational Mathematics, Department of Computational Mathematics, (e-mail: [email protected]). * Department of Computational Mathematics, Institute of of Computational Modeling Modeling of of SB RAS, RAS, Krasnoyarsk, Krasnoyarsk, Russia Russia Institute **Department of Mathematical SupportofforSB Devices andRussia Systems, Institute of Computational Computational SBDiscrete RAS, Krasnoyarsk, (e-mail: Modeling [email protected]). (e-mail: [email protected]). Siberian Federal University, Krasnoyarsk, Russia (e-mail: [email protected]). **Department of Mathematical Support for Discrete Devices and Systems, **Department of Support for Discrete and (e-mail: [email protected]) **Department of Mathematical Mathematical Support for Discrete Devices Devices and Systems, Systems, Siberian Federal University, Krasnoyarsk, Russia Siberian University, Krasnoyarsk, Russia Siberian Federal Federal University, Krasnoyarsk, Russia (e-mail: [email protected]) (e-mail: (e-mail: [email protected]) [email protected]) Abstract: An inequality for stability control for the Ceschino method of the second order is constructed. The First order method with an extended stability domain based on of stages of the order Ceschino scheme is Abstract: An inequality for stability control for the Ceschino method the second is constructed. Abstract: An inequality for stability control for the Ceschino method of the second order is constructed. constructed. Calculations results of stiff problems modeling, confirming the increased efficiency of Abstract: An inequality for stability control for the Ceschino method of the second order is constructed. The First order method with an extended stability domain based on stages of the Ceschino scheme is The First order method with an extended stability domain based on stages of the Ceschino scheme is numerical calculations due to using alternating order are given. The First order method with an extended stability domain based on stages of the Ceschino scheme is constructed. Calculations results of stiff problems modeling, confirming the increased efficiency of constructed. Calculations Calculations results results of of stiff stiff problems problems modeling, modeling, confirming confirming the the increased increased efficiency efficiency of of constructed. numerical calculations due to to using alternating alternating ordercontrol, are given. given. © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: explicit methods, and stability stiff problems. numerical due using order numerical calculations calculations due toaccuracy using alternating order are are given. Keywords: explicit methods, accuracy and stability control, stiff problems. Keywords: Keywords: explicit explicit methods, methods, accuracy accuracy and and stability stability control, control, stiff stiff problems. problems. Jacobi matrix and, consequently, to the increase of 1. INTRODUCTION computational second approach is based on the Jacobi matrix matrix costs. and, The consequently, to the the increase of Jacobi and, consequently, to increase of 1. INTRODUCTION Jacobi matrix and, consequently, to the increase of λ of the Jacobi estimation of the maximum eigenvalue The necessity of applying algorithms based on explicit 1. INTRODUCTION max computational costs. The second approach is based on the 1. INTRODUCTION computational costs. second approach is the computational costs. The The secondusing approach is based based on on methods appears inapplying numerical large-scale stiffonproblems matrix by the power iterations the λincrements of the the the Jacobi estimation of the the maximum eigenvalue The necessity of algorithms based explicit max of λ of the Jacobi estimation of maximum eigenvalue The necessity of applying algorithms based on explicit max λmax of the Jacobi estimation of the maximum eigenvalue solving as shown in [1–3]. Integration methods on a base of The necessity of applying algorithms based on explicit right-hand part of a system of ordinary differential equations methods appears in numerical large-scale stiff problems matrix by the power iterations using the increments of the methods in numerical large-scale stiff problems matrix the iterations using the increments of the implicit orappears semi-implicit numerical schemes usually use the methods appears in[1–3]. numerical large-scale stiff problems matrix bysubsequent the power power iterations using the increments of [7]. the with theby control of ordinary the inequality h ⋅ λmax equations ≤D solving as shown in Integration methods on a base of right-hand part of a system of differential solving as shown in [1–3]. Integration methods on a base of right-hand part of a system of ordinary differential equations decomposition of the Jacobi matrix [2], that is a separate solving as shown in [1–3]. Integration methods on a base of right-hand part of a system of ordinary differential equations estimation doescontrol not leadoftothe theinequality increase ofh computational implicit or semi-implicit numerical schemes usually use the This with the subsequent ⋅λ ≤ D implicit or schemes usually use with subsequent control of h D [7]. [7]. time-consuming task. In numerical suchmatrix a situation, the algorithms implicit or semi-implicit semi-implicit numerical schemes usually use the the costs max ≤ with the the subsequent control of the the inequality inequality h ⋅⋅ λλmax in all the considered situations [3–5, 7]. max ≤ D [7]. decomposition of the Jacobi [2], that is a separate decomposition of the Jacobi matrix [2], that is a separate This estimation does not lead to the increase of computational based on explicit formulas are more efficient if stiffness of decomposition oftask. the Jacobi matrix [2], that is algorithms a separate This estimation does to of time-consuming In such aa situation, This estimation does not not lead lead to the the increase increase of computational computational time-consuming task. In an such situation,tothe the algorithms costs in all the considered situations [3–5, 7]. the problem allows to get approximation a solution in a time-consuming task. In such a situation, the algorithms In this paper, the stability control inequality costs in all the considered situations [3–5, 7]. based on on explicit explicit formulas formulas are are more more efficient efficient if if stiffness stiffness of of costs in all the considered situations [3–5, 7]. for the second based reasonable time [3]. based on explicit formulas are more efficient stiffnessinofa order Ceschino method [8] is constructed using the approach the problem allows to get an approximation to aaifsolution In the stability control inequality for the the problem allows to get an approximation to solution in In this thisinpaper, paper, thefirst stability control inequality forstages the second second the problem allows to get an approximation to a solution in aa In given [7]. The order method based on the of the this paper, the stability inequality for second reasonable time [3]. order Ceschino method [8] control is constructed using thetheapproach approach The step control algorithm is usually based on accuracy reasonable time [3]. order Ceschino method [8] is constructed using the reasonable time [3]. numerical Ceschino formula with a stability interval extended order Ceschino method [8] is constructed using the approach in [7]. The first order method based on the stages of the control of control a numerical scheme. It is natural, because the given given [7]. The first order based the of to 32 in along the real axis constructed. Ainterval newstages integration The step algorithm is usually based on accuracy given in [7]. The firstformula orderismethod method based on on the stages of the the The step control algorithm is usually based on accuracy numerical Ceschino with a stability extended accuracy the approximation the main of numerical The stepofof control algorithm is is usually basedrequirement on accuracy Ceschino formula with aa stability interval extended algorithm of alternating order and step is formulated. The control a numerical scheme. It is natural, because the numerical Ceschino formula with stability interval extended control of numerical scheme. It is natural, because the to 32 along the real axis is constructed. A new integration calculations. However, applying algorithms control ofofaa the numerical scheme.isthe Itthe isintegration natural, because the to 32 the real is constructed. A new integration numerical results of axis stifforder problems solving confirming the accuracy approximation main requirement of to 32 along alongof the real axis is constructed. A formulated. new integration accuracy of the approximation is the main requirement of algorithm alternating and step is The based on explicit methods to stiff problems leads to the loss accuracy of the approximation is the main requirement of algorithm of alternating order and step is formulated. The increase of the efficiency due to alternating order and the calculations. However, applying the integration algorithms algorithm of alternating order and step is formulated. The calculations. However, applying the integration algorithms results of stiff problems solving confirming the of efficiency and reliability [4–5]. The leads reason is loss the numerical calculations. However, applying the integration algorithms numerical results of stiff problems solving confirming the stability control are given. based on explicit methods to stiff problems to the numerical results of stiff problems solving confirming the based on on explicit explicit methods to stiff stiff problems problems leadsonto to the loss increase of the efficiency due to alternating order and the contradiction between the accuracy and stability settling based methods to of efficiency due of efficiency efficiency and reliability [4–5]. The leads reasona the is loss the increase increase of the the are efficiency due to to alternating alternating order order and and the the of and reliability [4–5]. The reason is the stability control given. region which leads to a large amount of the recomputed of efficiency and reliability [4–5]. The reason is the stability control control are are given. contradiction between the accuracy and stability on aa settling stability 2. given. CESCHINO METHOD contradiction between the accuracy and stability on settling solutions, whereas thetothe step is much lessstability thanthethe maximum contradiction between and on a settling region which leads aa accuracy large amount of recomputed region which leads to large amount of the recomputed 2. CESCHINO CESCHINO METHOD allowable one. The additional stability control of a numerical region which leads to a large amount of the recomputed We consider the Cauchy problemMETHOD 2. solutions, whereas the step is much less than the maximum 2. CESCHINO METHOD solutions, whereas the step is much less than the maximum scheme allows to avoid this problem. At the present moment, solutions, whereas the step is much less than the maximum allowable one. The additional stability control of a numerical We consider the Cauchy problem y ′ = f (t, y ) , the y (t0 ) = y0 , problem t0 ≤ t ≤ t k , (1) allowable one. The additional stability of aa numerical there are two common approaches topresent stability We allowable one. The additional stabilityAtcontrol control ofthe numerical We consider consider the Cauchy Cauchy problem scheme allows to avoid this problem. the moment, scheme allows to avoid this problem. At the present moment, control [5–7], obtained applying the Dahlquist test equation. scheme allows to avoid this problem. At thetopresent moment, yy ′′ = f (ty, yand ) , y (ft ) are = yy0 N , t-dimensional ≤ tt ≤ t , (1) where real vector functions, there are two common approaches the stability (1) there are two common approaches to the stability y′ = = ff ((tt,, yy )) ,, yy ((tt000 )) = = y00 ,, tt000 ≤ ≤t ≤ ≤ ttkkk ,, (1) there are two common approaches to the stability control [5–7], obtained applying the Dahlquist test equation. t is an independent variable, which varies over the interval where and y f are N -dimensional real vector functions, The first approach refers with the estimating the maximum control [5–7], [5–7], obtained obtained applying applying the the Dahlquist Dahlquist test test equation. equation. where and y f are N -dimensional real vector functions, control and f are N -dimensional real vector functions, eigenvalue of the refers Jacobi with matrix its norm with the where For the solution of (1) we use theover explicit Runge[ttt0isis, tkan ] . yindependent variable, which varies the interval The first first approach approach the using estimating the maximum maximum an independent The refers with the estimating the t is an independent variable, variable, which which varies varies over over the the interval interval subsequent (in addition to the accuracy control) control of the The first approach refers with the estimating the maximum Kutta eigenvalue of the Jacobi matrix using norm with the . For the the solution of of (1) we we use the the explicit RungeRunge[[ttt00 ,,, tttkk ]]formulas eigenvalue of the matrix using its its with the eigenvalue of faddition the Jacobi Jacobi matrix its norm norm with the For the solution solution of (1) (1) we use use the explicit explicit Runge[ 0 k ] .. For ancontrol) integration stepof and inequality h isusing h(in subsequent to the accuracy control the y ≤ D , where subsequent (in addition to the accuracy control) control of the Kutta formulas yn +1 =formulas yn + pm1k1 + pm 2 k2 + pm 3 k3 + pm 4 k4 , subsequent (in addition to the accuracy control) control of the Kutta Kutta formulas ,, where is an integration step and inequality h h ff yof≤ D D is a norm ODE system Jacobi matrix [6]. Here, is f where is an integration step and inequality h h ≤ D y y inequality h f y ≤ D , where h is an integration step and y = yyn + pm1kk1 + pm 2 kk2⎛+ p 1k + pm 41kk4 ,, ⎞ 2 k2 + m 4 4k kyy1nnn +++=111 = hf y(nntn+ t p+mm 33 kh33 ,+y p+ , (2) = +, yp pnmm)11k, 11 k+ +2 p p=mmhf 2 2⎜+n pm 3 k3 + npm 4 k 4 ,1 ⎟ a ffpositive constant correlated with the size[6]. of Here, a stability D is aa norm of ODE system Jacobi matrix is y D is norm of ODE system Jacobi matrix [6]. Here, is 4 4 ⎞⎠ ⎝ 1 1 y ⎛ D is a norm of ODE system Jacobi matrix [6]. Here, is f 1 1 y ⎛ ⎞ domain. It is = hf hf ( t , y ) , k =hf hf ⎛⎜ ttn + + 1 h,, yyn + + 1 kk1 ⎞⎟ ,, (2) kkk11 = (2) positive constant constant correlated correlated with with the the size size of of aa stability stability = hf ((⎛ttnn ,, yynn1)) ,, kk22 = =hf h, ynn + 4 1 ⎜ t⎞n + 4 h (2) aaawell-known, positive 4 k11 ⎟⎟⎠⎠ , that explicit numerical schemes involve positive constant correlated with the size do of not a stability k13 =hf ⎜ tnn +n h, y2n + k⎜⎝⎝⎝2 ⎟n , 4 4 4 domain. It is ⎠ domain. It is 2 2 ⎠⎞ the Jacobi that matrix. Therefore, thenotapproach domain. It applying is ⎝⎛ 1 well-known, explicit numerical schemes do involve 1 1 ⎛⎜ tn + 1 kk3 = hf ,, yyn + well-known, that explicit numerical schemes do not involve 1h 1 kk2 ⎞⎞⎟ ,, ⎛ = hf t + h + mentioned above leads to the additional calculation of the well-known, that explicit numerical schemes do not involve k33 =hf ⎜⎝⎜ tnn + 2 h, ynn + 2 k22 ⎟⎠⎟ , the Jacobi matrix. Therefore, applying the approach 2 2 the matrix. Therefore, applying the ⎝⎝ 2 2 ⎠⎠ the Jacobi Jacobiabove matrix. Therefore, applying the approach approach mentioned leads to the additional calculation of the mentioned above leads mentioned© 2015, aboveIFAC leads to to the the additional additional calculation calculation of of the the 198 Copyright 2405-8963 © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Peer review©under of International Federation of Automatic Control. Copyright 2015, responsibility IFAC 198 Copyright © 198 10.1016/j.ifacol.2015.05.080 Copyright © 2015, 2015, IFAC IFAC 198
MATHMOD 2015 February 18 - 20, 2015. Vienna, Austria Eugeny A. Novikov et al. / IFAC-PapersOnLine 48-1 (2015) 198–203
k= hf ( tn + h, yn + k1 − 2k2 + 2k3 ) , 4
where X = h ⋅ A . It is easy to see, that the following relations
where h is an integration step, ki , 1 ≤ i ≤ 4 , are the stages of the method, pmi , 1 ≤ i ≤ 4 , are the numerical coefficients, m is the accuracy order of the method. On the coefficients p21 = 1 , p22 = −2 , p23 = 2 , p24 = 0 ,
199
k1 − 2k2 + k3=
(3)
the scheme (2), (3) has the second accuracy order [8].
1 3 X ⋅ yn , 8
1 1 2 k1 ) X ⋅ yn , ( k2 −= 2 8
hold. The estimation of the maximum eigenvalue of the Jacobi matrix of the system (1) can be calculated by the power iterations [3]. We introduce the notation
3. ACCURACY CONTROL The scheme (2) has the fourth order on the coefficients 1 2 1 , p 42 =0 , p 43 = , p 44 = . The estimation of the local 6 3 6 error δn ,2 of the form p 41 =
⎪⎧ ( k1 − 2k2 + k3 )i ⎪⎫ vn = 2 ⋅ max ⎨ ⎬. 1≤ i ≤ N ⎪⎩ ( k2 − k1 )i ⎪⎭
δn ,2 = ( p41 − p21 ) k1 + ( p42 − p22 ) k2 +
For stability control of the Ceschino method the inequality vn ≤ D can be applied, where the value of D bounds the stability interval.
+ ( p43 − p23 ) k3 + ( p44 − p24 ) k4 can be used for the accuracy control. As a result, for the accuracy control we use the inequality
(4)
The stability properties of the Runge-Kutta methods are usually investigated on the scalar test equation y ′ = λy , where λ is some arbitrary complex number, Re ( λ ) < 0 . λ
δn ,2 ≤ ε ,
pretends to be an eigenvalue of the Jacobi matrix of the problem (1). Applying (2), (3) to the solution of y ′ = λy we
where ⋅ is some norm in R N , ε is the required accuracy of calculations. Here, the relation δn ,2 = O ( h 3 ) holds.
get that the stability function Q2 ( x ) of the second order method has the form
In this paper the term "accuracy step" is referred to the step chosen according to the accuracy requirements and the term "stability step" denotes the step, which is chosen according to the stability requirements. We choose the accuracy step hac using the formula h ac= q ⋅ h , where q is determined from
Q2 ( x ) =1 + x +
1 2 1 3 x + x , x= h ⋅ λ , 2 4
and the stability function Q4 ( x ) of the fourth order method
the equation q 3 ⋅ δn ,2 = ε . If q < 1 , we recompute the
has the following form
solution with the step h equal to q ⋅ h . Otherwise, we compute and accept the solution and the predicted step hn +1 is calculated using the formula hn +1= q ⋅ h . The inequality
Q4 ( x ) =1 + x +
1 2 1 3 1 4 x + x + x , x= h ⋅ λ . 2 6 24
The length of the stability interval of the second order method is equal to 2 (fig. 1), whereas for the fourth order scheme it is equal to 2.8 (fig. 2). Therefore, in the inequality vn ≤ D we let D = 2 . Taking into account, that vn = O ( h ) ,
δn ,2 ≤ ε is well-behaved in solving a lot of practical
problems, therefore, we use it in this paper. In this paper the algorithm of alternating step based on the scheme (2), (3) with the inequality for accuracy control δn ,2 ≤ ε will be
the stability step h st can be determined from the formula hst = r ⋅ h , where r is calculated from the equation , D 2. r= ⋅ vn D=
referred to as CESCH42 .
There are three reasons why estimation (4) is rather rough: 1) the maximum eigenvalue may be very close to other eigenvalues; 2) there are only a few iterations in the power iterations method; 3) as the estimation is obtained on the linear problem y ′ = Ay , non-linearity of the problem (1) leads to the additional error in the estimation of the maximum eigenvalue. Since the estimation is rough the stability control is used to limit stepsize. As a result, the predicted stepsize is calculated using the formula
4. STABILITY CONTROL of NUMERICAL SCHEME Now we construct the inequality for stability control of the numerical scheme (2). For this purpose, we apply (2) for the solution of the linear problem y ′ = Ay with some constant matrix A . The first three stages k1 , k2 and k3 applied to this problem have the following form
1 2⎞ ⎛ k= X ⋅ yn , k 2 = 1 ⎜ X + X ⎟ ⋅ yn , 4 ⎝ ⎠ 1 1 ⎛ ⎞ k3 =⎜ X + X 2 + X 3 ⎟ ⋅ yn , 2 8 ⎝ ⎠
{
}
hn +1 = max hn , min ⎡⎣ h ac , h st ⎤⎦ ,
199
(5)
MATHMOD 2015 February 18 - 20, 2015. Vienna, Austria Eugeny A. Novikov et al. / IFAC-PapersOnLine 48-1 (2015) 198–203 200
settling region is much higher than the required one. It is natural, because on this region the previous errors are suppressed due to the stability control, whereas new errors are small due to low values of the derivatives of the solution. In this situation, it is more efficient to perform calculations using some low order method with wider stability domain.
where hn is the latest successful integration step. Notice, that formula (5) is applied to predict the stepsize of hn +1 after successful calculation of a solution with the previous step hn and, therefore, does not lead to the increase of computational costs. If the stability step is less than the latest successful one, it will not be decreased, because the demand of decreasing the stepsize may be caused by rather roughness of the estimation. However, the stepsize will also not be increased, because the possibility of instability of a numerical scheme is not excluded. If the stability step has to be decreased, the latest successful step hn is applied. As a result, the formula (5) is used to choose the stepsize. The given formula allows to stabilize the step behavior on the settling region of the solution, where the stability plays a defining role. It is the existence of the settling area, that significantly restricts the capabilities of the explicit methods on solving stiff problems.
5. FIRST ORDER METHOD Now we construct the first order method yn +1,1 = yn + p11k1 + p12 k 2 + p13 k3 + p14 k 4 ,
(6)
based on the stages of the numerical scheme (2) with a wider stability domain. For this purpose, we apply (6) to the solution of the test equation y ′ = λy . We get
yn +1,1 = Q1 ( x ) yn , where x = h ⋅ λ and the stability function Q1 ( x ) has the form Q1 ( x ) = 1 + ( p11 + p12 + p13 + p14 ) x + 1 ⎛1 ⎞ + ⎜ p12 + p13 + p14 ⎟ x 2 + 2 ⎝4 ⎠ 1 1 ⎛1 ⎞ + ⎜ p13 + p14 ⎟ x3 + p14 x 4 . 8 2 4 ⎝ ⎠
(7)
The requirement of the first accuracy order of the numerical formula (7) leads to the relation p11 + p12 + p13 + p14 = 1.
We use the other coefficients p1i , 1 ≤ i ≤ 4 , to make the stability domain wider. The stability condition for the scheme (6) has the form Q1 ( x ) ≤ 1 . In order to construct a method
Fig. 1. The stability domain of the second order method.
with the maximum stability interval we consider the Chebyshev polynomial
T4 ( z ) = 8 z 4 − 8 z 2 + 1 . It is known, that the polynom T4 ( z ) is less deflected from zero for z ∈ [ −1,1] . We introduce the notation z= 1−
2 x, γ
so that the interval [ γ, 0] is mapped to the interval [ −1,1] . As a result, the polynom T4 ( x ) has the form T4 ( x ) = 1−
32 160 256 128 x + 2 x 2 − 3 x3 + 4 x 4 . γ γ γ γ
It is easy to show [3], that for the polynom T4 ( x )
Fig. 2. The stability domain of the fourth order method.
(8) the
inequality T4 ( x ) ≤ 1 is satisfied in the maximum interval
In this paper, the algorithm of alternating step with the additional stability control of the numerical scheme will be referred to as CESCH42st . The given algorithm is based on a numerical formula of low (the second) accuracy order, therefore, it is aimed at the solution of non-stiff and moderate stiff problems with low accuracy of calculations (of 1% order and lower). From the calculations with the algorithm CESCH42st it follows, that the practical accuracy on a
[ γ, 0 ] ,
γ = −32 (fig. 3). Comparing the relations (7) and (8)
on γ = −32 , we get the coefficients p1i , 1 ≤ i ≤ 4 , of the first order method (2) with the maximum stability interval, i.e. p11 = 200
895 257 31 1 , p12 = , p13 = , p14 = . 2048 512 512 2048
(9)
MATHMOD 2015 February 18 - 20, 2015. Vienna, Austria Eugeny A. Novikov et al. / IFAC-PapersOnLine 48-1 (2015) 198–203
In the fig. 3 the stability domain consists of four parts, which intersect in extreme points of the Chebyshev polynomial. For those problems, which Jacobi matrix eigenvalues have nonzero imaginary part, it may lead to decrease of the calculation efficiency. This effect may be overcome by applying stability polynomials which are close to optimal polynomials. In [3] there is given an algorithm for construction stability polynomials with stability domain extended along the implicit axis.
201
order method are accompanied by (in addition to accuracy control) the control of the inequality vn ≤ 32 , whereas the stepsize is chosen using the formula (5). In this paper the algorithm of alternating order and step is referred to as CESCH42vp . 8. DIFFERENTIAL EQUATIONS of CHEMICAL KINETICS The kinetic scheme of any chemical reaction involves the following elementary reactions (stages) α1,1 x1 + ... + α NR,1 xNR → β1,1 x1 + ... + βNR,1 xNR ,
……… α1,NS x1 + ... + α NR,NS xNR → β1,NS x1 + ... + βNR,NS xNR ,
(10)
where xi , 1 ≤ i ≤ NR are chemical reagents, NR and NS are the number of the reagents and the stages in the reaction, respectively; αij and βij , 1 ≤ i ≤ NR , 1 ≤ j ≤ NS are the stoichiometric coefficients. For each elementary reaction the corresponding stage velocity constant k j , 1 ≤ j ≤ NS is
Fig. 3. The stability domain of the first order method.
given. For the process (10) under the lumped model of the constant-volume isothermal reactor the corresponding ordinary differential equations system has the following form
The stability domain of the first order method (6), (9) along the real axis is 16 times wider than the stability domain of the numerical scheme (2), (3). The computational costs of the first and the second order method are the same. Therefore, for the problems, in which the stepsize is limited mostly due to the stability issues the assumed theoretic increase of the efficiency is equal to 16 times.
C ′ = AT V , C ( 0 ) = C0 . Here, AT is a stoichiometric matrix, C and V are vectors of reagents concentrations and stage velocities, respectively. In case of the reaction flowing in non-isothermal conditions this system also involves the heat balance equation
6. ACCURACY AND STABILITY CONTROL of FIRST ORDER METHOD
T′ =
In the accuracy control inequality we use the following estimation of the local truncation error
QT V − α (T − T01 ) CVT C
,
where T is the mixture temperature in the reactor, T01 is the
δn ,1 = ( p41 − p11 ) k1 + ( p42 − p12 ) k2 +
temperature of the reactor walls, QT is the vector of the
+ ( p43 − p13 ) k3 + ( p44 − p14 ) k4 .
calorific capacities of the stages, CVT is a vector of the heat % / r , α% is the heat-conduction capacities of the reagents, α = αs coefficient, s and r are square of area and reactor volume, the vectors respectively. The superscript T for QT and CVT denotes the transposition. The heat capacities of the reagents and the heat-conduction coefficient may be functions of the reagents concentrations ci , 1 ≤ i ≤ NR , α may depend on the temperature.
For the accuracy control of the numerical formula (6), (9) we can apply the inequality δn ,1 ≤ ε . Since the stability interval of the numerical scheme (6), (9) is limited by 32, we apply for stability control the inequality vn ≤ 32 where vn is determined using the formula (4). 7. ALGORITHM of ALTERNATING ORDER
If the reaction proceeds in a constant-volume isothermal reactor with the substance exchange (an open system, a perfect-mixing reactor) a system of ordinary differential equations has the form
The first order method with the extended stability domain is efficient on the settling regions, where stepsize is limited due to the stability issues. On the transition region it is more efficient to apply the method (2), (3). We can increase the efficiency by applying each method on the region, where it is the most efficient. As a criteria for switching between the methods we can use the inequality for stability control. On calculations with the method (2), (3) we switch to the numerical scheme (6), (9) on the violation of the inequality vn ≤ 2 . On calculations with the first order method we switch back when vn ≤ 2 satisfies. The calculations by the first
C ′ =AT V +
1 ( C p − C ) , C ( 0 ) = C0 . Θ
where C p is a vector of incoming reagent concentrations, Θ is the stay period of the mixture in the reactor, Θ = r / u , u is the mixture space velocity through a reactor. If the
201
MATHMOD 2015 February 18 - 20, 2015. Vienna, Austria Eugeny A. Novikov et al. / IFAC-PapersOnLine 48-1 (2015) 198–203 202
10. NUMERICAL MODELING OF ETHANE PYROLYSIS
reaction proceeds in non-isothermal conditions this system also involves the heat balance equation
T′
QT V − α (T − T01 ) CVT C
The ethane pyrolysis in absence of oxygen is determined by a small sequence of elementary reactions. The mechanism of the ethane pyrolysis has been discussed in scientific literature many times. Here, we use the reaction scheme, which was proposed and investigated in [10]
1 − (T − T02 ) , Θ
where T02 is the incoming mixture temperature in the reactor. The temperature of the reaction mixture can be described by a function of time t and concentrations of reagents ci ,
C2 H 6 → CH 3 +CH3 , CH 3 +C2 H 6 → CH 4 +C 2 H 5 , C2 H 5 → C2 H 4 +H , H+C 2 H 6 → H 2 +C2 H 5 ,
1 ≤ i ≤ NR , i.e. T = T ( t , C ) .
C 2 H 5 +C2 H 5 → C 4 H10 .
9. ALGORITHM FOR GENERATING CHEMICAL KINETIC EQUATIONS
Here, the stage velocity constants have the form = k1 1.34 ⋅ 10−5 , k 2 = 3.73 ⋅10 2 ,= k3 3.69 ⋅103 ,
If an elementary reaction is balancing, its velocity Ws is
k 4 = 3.66 ⋅ 105 и k5 = 1.62 ⋅107 .
equal to a difference of its straight Ws+ and backward Ws−
We denote reagent concentrations in the following way
reactions, i.e. W = Ws+ − Ws− , 1 ≤ s ≤ NS . If there is some s third particle in a reaction, the velocity Vs is calculated using the formulas [9]
c1 = [ C2 H 6 ] , c2 = [ CH3 ] , c3 = [ CH 4 ] , c4 = [ C2 H5 ] ,
Vs = PW s s , Ps =
NR+NI
∑ i =1
c5 = [ C2 H 4 ] , c6 = [ H ] , c7 = [ H 2 ] , c8 = [ C4 H10 ] . Then, the corresponding system involves 8 ordinary differential equations of the form C ′ = AT V , i.e.
εsi ci , 1 ≤ s ≤ NS ,
c1′ = −k1c1 − k2 c1c2 − k4 c1c6 = , c2′ 2k1c1 − k2 c1c2 ,
where εsi , 1 ≤ s ≤ NS , NR + 1 ≤ i ≤ NR+NI , are the efficiencies of the third particles, NI is the number of inert substances, εsi and ci are the efficiencies and concentrations of the inert substances, respectively. The values of components of the vector Ws are determined from the scheme of chemical reaction (10) using the relations + s
W = ks
NR+NI
Πc i =1
αij
− s
, W = k− s
NR+NI
Πc i =1
βij
c3′ = k2 c1c2 , c4′ = k 2 c1c2 − k3 c4 + k 4 c1c6 − 2k5 c42 ,
(11)
c5′ = k3 c4 , c6′ = k3c4 − k4 c1c6 , c7′ = k4 c1c6 , c8′ = k c . 2 5 4
The initial concentration of the ethane c1 = [ C2 H 6 ] is equal to 0.14, other reagent concentrations are equal to zero. The calculations were performed with the accuracy ε = 10−2 . The efficiency of the integration algorithms was estimated by the total number if of calculations of the right-hand part of the problem (11) over an integration interval. The numerical approximation was performed over the interval [0,0.26]. The given problem satisfies the 'classical' definition of stiffness. In the beginning of the interval a transition region (centiseconds) is observed and after there is a settling region (fig. 4).
,
where k s and k− s , 1 ≤ s ≤ NS , are straight and backward velocity stage constants, respectively. Velocity constants of the stages are calculated using the formulas ⎛ Ej ⎞ n = k j Aj T j exp ⎜ − ⎟, ⎝ RT ⎠ where T is the temperature of the mixture in a reactor; Aj , n j and E j / R are the given constants. Notice that, in
general, the velocity constants are not constant in case of a non-isothermal reactor – they depend on temperatures. However, historically, an isothermal reactor was considered earlier than the non-isothermal one and thus k j , 1 ≤ j ≤ NS , at the present moment, are still called constants. The stoichiometric matrix AT with the elements aij is formed from the kinetic scheme (10) conforming to the following rule. The number of a stage is aligned with the column number and the reagent number is aligned with the row number of the matrix AT . If xi serves as an initial reagent, then aij = − αij , if xi is a product, then aij = βij . If xi is both an initial reagent and a product, then aij = −αij + βij . Usually a few amount of reagents are reacting in a stage, i.e. a stoichiometric matrix is sparse.
Fig. 4. Time-dependence of CH 3 concentration.
202
MATHMOD 2015 February 18 - 20, 2015. Vienna, Austria Eugeny A. Novikov et al. / IFAC-PapersOnLine 48-1 (2015) 198–203
The comparism of the efficiencies was performed between the constructed algorithms and the known Merson method (MERSON, [11]). In case of applying explicit methods for stiff problems solving the integration stepsize is limited due to the stability requirements almost over the whole integration interval. As a result, the practical accuracy of a solution is much higher than the defined one. In fig. 4 the graphs of a solution are coincide for the Merson method and for the constructed algorithm CESCH42.
203
6. Shampine, L.M. (1982). Implementation of Rosenbrock methods. ACM Transaction on Mathematical Software, No. 5, 93–113. 7. Novikov, V.A. and Novikov, E.A. (1984) Control of the stability of explicit one - step methods of integrating ordinary differential equations. Soviet Math. Dokl., vol. 30, No. 1, 211–215. 8. Ceschino, F., and Kuntzman, G. (1966) Numerical solution of initial value problems. Prentice-Hall, Englewood Clis, New Jersey. 9. Novikov, E.A. (2010) Numerical modeling of a modified oregonator by the (2,1)-method for solving stiff problems. Numerical methods and programming, vol. 11, No. 2, 123–130 (in Russian). 10. Kulich, D.M., and Taylor, J.E. (1975) Mathematical simulation of the oxygen ethane reaction. Chem. Kinetic, vol. 8, 89–97. 11. Merson, R.H. (1957) An operational methods for integration processes. Proc. of Symp. on Data Processing, Australia, 329–331. 12. Enright, W.H., and Hull, T.E. (1975) Comparing numerical methods for the solutions of systems of ODE’s. BIT, No. 15, 10–48.
For CESCH42 without the stability control we have if = 22 853 , for the algorithm CESCH42st with stability control we have if =20 403 , for the algorithm of alternating order and step CESCH42vp the amount of calculations of the right-hand part is if = 2 588 and for the Merson method it is equal to if = 25 796 . 11. CONCLUSIONS From the numerical results we may draw the following important conclusions. Firstly, the constructed integration algorithm of the second order with accuracy control of calculations and the stability control of the numerical scheme as well as the algorithm of alternating order and step can be applied for the solution of rather stiff problems. Secondly, by the calculation costs the algorithm of alternating order and step CESCH42vp is more efficient than the Merson method almost 10 times, that is due to the stability control of the numerical scheme and calculations with alternating order. It seems that at rather high dimension of the problem (1) the method CESCH42vp may compete with implicit methods at solving moderate stiff problems, because it does not require to decompose the Jacobi matrix. The algorithm CESCH42vp is more efficient at the solution of the twelve test problems [2] and ten examples [12]. ACKNOWLEDGMENT This work was partially supported by Project 14-11-00147 of Russian Scientific Foundation. REFERENCES 1. Hairer, E., Norsett, S.P., and Wanner, G. (1987). Solving ordinary differential equations. Stiff and differentialalgebraic problems, Springer-Verlag, Berlin. 2. Hairer, E., and Wanner, G. (1996). Solving ordinary differential equations. Non stiff problems, SpringerVerlag, Berlin. 3. Novikov, E.A. (1997). Explicit methods for stiff systems, Nauka, Novosibirsk (in Russian). 4. Novikov, E.A., and Shornikov, Yu.V. (2013). Computer simulation of stiff hybrid systems,: NSTU Publishers, Novosibirsk (in Russian). 5. Novikov, E.A., and Zacharov, A.A. (2011). Harmonization of the stability regions in the three-stage explicit RungeKutta. Vestnik TyumGU, No. 7, 187–192 (in Russian).
203