Applied Mathematics
On a global error estimate
in long-term numerical integration
of ordinary dierential equations
G.A. Chumakov
1, N.A. Chumakova
2?1 Sobolev Institute of Mathematics, SB RAS, Pr. Ak.Koptyuga 4, Novosibirsk,
Russia, e-mail:chumakov@math.nsc.ru
2 Boreskov Institute of Catalysis, SB RAS, Pr. Ak.Lavrentieva 5, Novosibirsk,
Russia, e-mail:chum@catalysis.nsk.su
Received: March 06, 2001
Summary.
We describe an eective computational procedure of obtaining estimate for the global error in the long-term numerical integration with one-step embedded explicit Runge-Kutta methods applied to an ODE system. The paper relies on the existence of an asymptotic expansion for the global error. An example is given using a particular ODE system modeling the behavior of a heterogeneous catalytic reaction with complex dynamics.Key words:
numerical solution of ODEs, Runge-Kutta methods, global error, chemical kinetic modelMathematics Subject Classication (1991): 65-01, 65L06, 65Y99
1. Introduction
Mathematical models constructed in science and engineering are of-ten nonlinear and formed by complicated systems of dierential equa-tions. To study the dynamic behavior of such models one needs both the modern mathematical theory of dynamical systems and the prac-tical computational procedures (see some examples in 7{9]).
? This research was supported in part by the International Association for the
promotion of co-operation with scientists from the New Independent States of the former Soviet Union (INTAS) grant No. 99-01882.
In this paper we use the asymptotic expansion theory to analyze the global discretization error of numerical integration of the initial value problem
(1) z_=f(z) z(0) =z0 0 t T
where z and f(z) are m-dimensional vector functions and z0 is a constant vector. We assume that all components off are suciently smooth, in the following examples we assumefi (i= 1:::m) to be analytical.
Having dened tn=nh forn= 012:::and solved some nite-dierence equation that approximates (1), we would obtain an ap-proximation fx
n
g of the true solution z(t) at each t
n. We shall as-sume that the sequencefx
n
g the solution to the discrete problem -can be obtained exactly.
Note that any numerical integration is subject to roundo errors, but we assume that they are negligible comparingto the discretization error. Dening
en=xn ;z(t
n)
we shall call en the global discretization error at the end-point of the interval 0tn]. The aim of this paper is to produce some sort of estimate for en.
The term "error estimate" can have two meanings: either the error bound, or the error behavior description. The error boundary, though being the only way to obtain sharp estimates, must take into account the worst possible cases thus providing an often too pessimistic of a prediction. On the other hand, the sharp error estimate that may be several orders of magnitude bigger than the real error, is practically useless.
The subject of our investigation is to construct an asymptotic approximation of the function en as some small, e.g., h, is rened. By the term asymptotic approximation by the parameter h for the function en we mean the function e(th) such that the dierence
en ;e(t
nh), called "the remainder of the asymptotic approximation", is small for 0 tn T when the parameter h is small enough. The practical value of this approach is that we can eectively calculate
e(th) using simple dierential equations.
It should be noted that regardless of whether we want the precise sharp error estimate or the asymptotic approximation of the global error would suce, some auxiliary computations are necessary.
To analyze the global discretization error we consider the local errors of approximation as a result of perturbation of the original
initial problem (1). We regard the approximate solution x(t), where
x(tn) =xnas the exact solution of the initial value problem (2) xx_(0) ==fz(x) +(tx)
0
with the right hand side that diers from the initial problem (1) by the quantity(tx). Using the function(tx) we simulate the nite-dierence scheme, where(tx) is the local error of the numer-ical integration process, and is the small parameter that depends of the grid steph.
As the approximate value of (x;z) we will use the rst and the only term of the asymptotical decomposition
(3) x(t);z(t) = (t) + (t)
where (t) !0 as !0 and t 2 0T]. The function (t) is the solution of the initial value problem
(4) _ =f0
x(x(t))+(tx)
(0) = 0
where f0
x(x(t)) is the Jacobian matrix of rst partial derivatives of the functionf.
While calculating the local discretization error(tx) in (2) for the function(tx) we will also use the asymptoticexpansion ~(tx) such that the dierence(tx);~(tx) will be small for suciently small. In this case we can give an a priori estimate of;~, where ~
is the solution of the initial value problem (5) ( _~ =f0 x(x(t)) ~+ ~(tx) ~ (0) = 0 Let f0 x(x(t)) and
;~ be bounded on the interval 0 t T, namely
k(tx);~(tx)k N() kf 0 x(x(t))
k M
whereN()!0 when!0. Here the matrix norm is subordinated to the preassigned vector norm. Then;~on the given time interval satises the inequality
(6) k;~k k(;~)(0)ke Mt+N() M (eMt ;1) = N() M (eMt ;1):
It follows from (6) that k;~k !0, when ! 0. It should be noted that in sti problem 1] the classical Lipschitz constant
M = sup 0tT kf 0 x(x(t)) k 1
and consequently using the inequality (6) to estimate the norm of the
;~on the interval 0T] where T 1 can lead to the unacceptably severe restrictions on the parameter.
To estimate numerically the global discretization error we will simultaneously integrate the systems (1) and (4). With that method we checked in the Sect. 5 the reliability of the numerical integration of a particular nonlinear dynamic system that models a heterogeneous catalytic reaction behavior 4{6].
2. Runge-Kutta method structure
Numerical integration of the initial value problem (1) consists of a sequential calculation of values xn+1 as a discrete problem solution at the pointstn+1 forn= 01:::
For numerical solution of the initial value problem (1) and (4) we will use R-K explicit methods and follow 1{3]. These methods are described by (7) xn+1=xn+h s X i=1 biKi where (8) Ki=f(tn+cihxn+h s X j=1 aijKj) 1 i s aij = 0 for j i: Sometimes we will denote P
s i=1b
iKi as (hx(tn)). The step h can change withn.
Formulas (7) and (8) are called s-stage R-K formula because it is based onscalculations of the function f. Coecients in formula (8) usually satisfy the conditions
(9) ci = s X j=1 aij 1 = s X i=1 bi:
Table 1.
s 1 2 3 4 5 6 7 8 9 p 1 2 3 4 4 5 6 6 7
Compact form of a R-K method is given in a so-called the Butcher table (10) c AbT = c1 a11::: a1s c2 a21::: a2s ... ... ... cs as1::: ass b1 ::: bs
Note that there exist explicit R-K methods of an arbitrary order 3]. However, there is no general result about the maximumobtainable order p for given number of stages,s. Explicit R-K s-stage methods of the order p = s exist only for s= 123 and 4, but there are no 5-th order methods for s=5 3]. Table 1 shows maximum order p of an s-stage explicit R-K method.
2.1. Local error
Letx(t) be the exact solution of the system _x=f(x) on the segment
tn t tn+1with the initial datax(tn) =xn. Let us dene the local discretization error on tntn+1] as the dierence between xn+1 and the exact solutionx(tn+1):
(11) n+1=xn+1
;x(t n+1):
Thus, if R-K method has the approximation order p then the local discretization error of the R-K method is O(hp+1) withh
!0.
2.2. Local error estimate
For the numerical computation we distribute the nodes to make nu-merical integration error roughly the same on every interval. For this we need to know the computational error on every integration step. Now we describe some strategies of its numerical estimate.
While solving the initial value problem
with the R-K method of orderp, the main term of the local discretiza-tion error is (13) '(p+1)(t n) h p+1 (p+ 1)! where the function '(p+1) does not depend onh.
Let h = (tn+1 ;t
n)=2. Points (tn+hx(tn+h)) and (tnxn) are close, consequently the error on the next step will have the same main term. In two steps we will obtain approximation xn+1 of x(tn+ 2h), such that (14) x n+1 ;x(t n+ 2h) 2 '(p+1)(t n) h p+1 (p+ 1)! :
If we apply the same R-K method with a step 2hat the point (tnxn), we will obtain the approximate value ofxn+1 for which
(15) xn+1 ;x(t n+ 2h) '(p+1)(t n) (2h) p+1 (p+ 1)! :
It follows from (14) and (15) that the main term of the local dis-cretization error for the step = 2h is
(16) n+1=xn+1 ;x(t n+ ) 2p 2p ;1 ( xn+1 ;x n+1): If, after the computation of xn+1
;x
n+1 we see that the local error is still too large, we should decrease the time step.
In order to choose the time step in a more intelligent way, we should be able to estimate the local discretization error while com-puting the least possible number of right hand sides f(x).
Another and more ecient way of calculating the local error con-sists of using a R-K method with (s+ 1) stages and order (p+ 1), with the Butcher table
(17) c1 a11::: a1s 0 ... ... ... ... cs as1::: ass 0 cs+1 b1 ::: bs 0 ~b1 ::: ~bs ~bs+1 aij = 0 for j i such that thes-stage method
(18)
c1 a11::: a1s ... ... ...
csas1::: ass
is ofp-order and is embedded in (17).
If we use the embedded method (18) to nd xn+1
(19) xn+1 =xn+h s X j=1 bj f(Yj) where (20) Yi =xn+h s X j=1 aij f(Yj) i= 1:::s and use the (s+ 1)-stage method (17) to calculate xn+1
(21) x n+1 =xn+h s+1 X j=1 ~ bj f(Yj) where (22) Ys+1=xn+h s X j=1 bj f(Yj)
then the estimate of the local discretization error of the method (18) will be (23) xn+1 ;x n+1=h s X j=1 (bj ;~b j) f(Yj) ;h~b s+1 f(Ys+1): In numerical examples in the Sect. 5 we used a 4-stage R-K method of the 3rd order for computation of the approximate solution. The method had the following Butcher table
(24) 0 0 0 0 0 1=3 1=3 0 0 0 1=3 1=6 1=6 0 0 1=2 1=8 0 3=8 0 1=2 0 ;3=2 2
To estimate the local discretization error we used the 5-stage R-K-M method of the 4th order
(25) 0 0 0 0 0 0 1=3 1=3 0 0 0 0 1=3 1=6 1=6 0 0 0 1=2 1=8 0 3=8 0 0 1 1=2 0 ;3=2 2 0 1=6 0 0 2=3 1=6
for which method (24) is embedded.
The order of a particular method can be found using Taylor de-composition ofxn+1
;x(t
n+1) where
(26) xn+1 =x(tn) +h (hx(tn)):
The equation (26) represents an one-step explicit R-K method for solution of the initial value problem (1).
To illustrate the statement about the order of approximation of methods (24) and (25), let us consider the case when (1) is a scalar equation. Using Taylor formula for the exact solution x(t) of initial value problem (12) we obtain
(27) x(tn+1) ;x n h =f+ 2!hff0+h 2 3! h f2f00+f(f0)2 i + +h3 4! h f3f000+ 4f2f0f00+f(f0)3 i + +h4 5! h f4f(4)+ 4f3(f00)2+ 7f3f0f000+ 11f2(f0)2f00+f(f0)4 i +O(h5) wheref and its derivatives are calculated at the pointxn.
For the approximate solution xn+1 calculated with (24), we have
(28) xn+1 ;x n h = 12K1 ; 3 2K3+ 2K4= =f+2!hff0+h 2 3! h f2f00+f(f0)2 i +h3 4! 7 9f3f000+ 3f2f0f00+f(f0)3 + +h4 5! 115 2333f 4f(4)+ 35 23 2f 3(f00)2+ 385 2233f 3f0f000+ 495 2332f 2(f0)2f00 + +O(h5)
For the approximate solutionxn+1calculated with (25), we have
(29) xn+1 ;x n h = 16 K1+ 4K4+K5] = =f+2!hff0+h 2 3! h f2f00+f(f0)2 i +h3 4! h f3f000+ 4f2f0f00+f(f0)3 i + +h4 5! 25 24f4f(4) + 256f3(f00)2+ 85 12f3f0f000 + 858f2(f0)2f00 + 56f(f0)4 +O(h5):
The coecientsKi(i= 1345) are determined by (8), (9) and (24), (25). From (27) and (28) it follows that for the 4-stage method (24)
n+1= xn+1 ;x(t
n+1) =O(h
p+1) p= 3
i.e., the method (24) is of the third order of approximation. From (27) and (29) it follows that for the 5-stage method (25)
n+1=xn+1 ;x(t
n+1) =O(h
p+1) p= 4 i.e., the method (25) is of the fourth order of approximation.
Since the method (24) is embedded in (25), we may consider the following as an estimate of the local discretization error of the method (24) and as the function (tx) in the system (2) on tntn+1]: (30) x n+1 ;x n+1 = h 3 K1 ; 9 2K3+ 4K4 ; 1 2K5 : 2.3. Grids
The choice of the nite-dierence grid plays the outmost important role in the asymptotic analysis of the discretization error of a R-K method. Let us concentrate on the choice of the gridfD
N g.
Consider the interval 0T] and the nodestn,n= 01:::N, such that
0 =t0 < t1 < ::: < tN;1 < tN =T:
We will call this set of pointsa grid and denote it asDN. Quantities
hn=tn ;t
n;1 (n= 1:::N) are calledthe grid steps of DN. The grid sequence fD
N g N2N 0 is called uniform if DN = ft n= T n N n= 01:::Ng N 2N 0 i.e., grid steps hn = T=N (n = 01:::N
;1) for every D N are constant.
For a non-uniform grid sequence fD N
g N2N
0, we can consider the following denition of the grid steps. The grid sequence fD
N g
N2N 0 is calledlinearly convergent if there existc1>0 andc2>0 such that every step hn in every gridDN satises the following:
c1 N hn c2 N = ^hn n= 01:::N ;1 N 2N 0: Then when N ! 1, quantities O(^h
p
n) and O(h p
n) are equivalent to
We shall choose on 0T] a linearly convergent sequencefD N
g N2N
0 such that there exists a piecewise constant function
: 0T] ! c 1c2] 0< c1 c2 such that hn = T N (tn;1) n= 1:::N:
Such is calledthe step length function for the grid sequence fD N g, and such fD N g N2N 0 is calledcoherent.
Example 3]
.Let the grid DN be an element of the sequence fD N g N2N 0, N 0 = f3m m2Ng. Let the nods ofD
N be tn = 3 =4T=Nn n= 0:::2=3N 3=2T=Nn;T=2 n= 2=3N+ 1:::N: Then hn= T N (t n;1) n= 1:::N where (t) = 3=4 t2 0T=2) 3=2 t2 T=2T]:
We will presume thatN denes theDN grid uniquely, the number of nodes inDN will be N+1 and considerN
;1 or a function ofT=N as a small parameterin system (2).
Usually in practical applications the grid DN together with the discretization parameterN are dened during the process of solution of the initial value problem. The recurrent denition ofDN is de facto a part of the discretization method.
2.4. The initial grid
To construct the initial grid let us nd out how small the time steps of the discretization problem must be. We will use the method stability function. It is dened by behavior of the 1-step method in the model problem
(31) y_=y 2C t0 y(0) =y 0:
In spite of the simplicity of the equation (31), it is broadly used as a model to predict the computational method stability. Any R-K method (7)-(8) applied to (31) can be written as
(32) yn+1 =R(z) yn z =h
whereR: C!C is a polynomial in case of the explicit methods, or a rational function with real coecients in case of the implicit R-K methods. The functionRis calledthe method stability function.
Set of pointsz2C such that
(33) jR(z)j 1
is calledthe area of the method's absolute stability.
If an absolute stability area includes the whole left half-plane Re(z) 0 then the method is calledA-stable. Unfortunately, the nec-essary condition of theA-stability is the implicity of the R-K method 3].
The explicit R-K method with the Butcher table (24) has the stability function in a polynomial form
(34) R(z) = 1 +z+ 12!z2
+ 13!z3
+ 14!z4:
There is a convenient formula 3] for calculation the stability function using the Butcher table (10)
(35) R(z) = det(I;zA;zeb T) det(I;zA)
whereeT = (1:::1).
Note that the stability functionR(z) for any R-K method of order
p satises the condition
(36) R(h) =eh+O(hp+1) for
jhj!0:
However (36) is a necessary, but not sucient condition for a R-K method to have the orderp. For example, considering (28) it follows that for the method (24) the stability function is
R(h) =eh+O(h5)
though the method (24) is of the third order of approximation. Consequently, while choosing the initial grid the time steps must be small enough to satisfy the condition
For example, it is sucient to require that hjj 1 or in case of the system (1) this condition can be written as
hjjf 0 x(x
n) jj 1:
Furthermore, maximum of the local error vector must be main-tained below the certain level (or in an"-neighborhood of it) on every step of the initial grid construction.
3. Choice of the function
In this section we describe how the small parameterand the func-tionof the system (4) are to be chosen when we consider the initial value problem (1).
Consider the rst time interval 0t1] on the initial grid DN. Let the main integration formula be of the order approximation of p. Thus if the solution z(t) of (1) on 0t1] has the form
(37) z( ) = p+1 X n=0 a0 n n+O( p+2) then the approximate solutionx(t) can be written as (38) x( ) = p X n=0 a0 n n+b0 p+1 p+1+O( p+2) moreoverb0 p+1 6 =a0 p+1 in general. From (37) and (38) it follows that (39) x( ) =z( ) + (b0 p+1 ;a 0 p+1) p+1+O( p+2):
Assumef(z) is smooth enough, and dierentiate (39) with respect to
. Taking into account that _z=f(z), we obtain (40) dxd =f x( );(b 0 p+1 ;a 0 p+1) p+1 ;O( p+2) + +(p+ 1)(b0 p+1 ;a 0 p+1) p+O( p+1): Decomposingf into series
0= (b 0 p+1 ;a 0 p+1) p+1+O( p+2)
and ignoring terms of orderO( p+1), we conclude that on the interval 0t1]x( ) satises (41) dxd =f(x) + (p+ 1)(b0 p+1 ;a 0 p+1) p:
Let us show that on the intervals titi+1], i = 1:::N
;1 of the gridDN, the function x(t) satises the similar equation. Let the solutionz(t) on the interval t1t2] be
(42) z(t1+ ) = p+1 X n=0 a1 n n+O( p+2) where =(t2 ;t 1) 2 01]: Then we can present the approximate solutionx(t) as
(43) x(t1+ ) = p X n=0 a1 n n+ +b1 p+1 p+1+ (b0 p+1 ;a 0 p+1)h p+1 1 +O(h p+2 1 ) +O( p+2): It follows from (42) and (43) that
(44) x(t1+ ) =z(t1+ )+ +(b1 p+1 ;a 1 p+1) p+1+ (b0 p+1 ;a 0 p+1)h p+1 1 +O(h p+2 1 ) +O( p+2): Dierentiating (44) with respect to and taking into account _z =
f(z), we get (45) dxd =f x(t1+ ) ;(b 1 p+1 ;a 1 p+1) p+1 ; ;(b 0 p+1 ;a 0 p+1)h p+1 1 +O(h p+2 1 ) +O( p+2) + +(p+ 1) (b1 p+1 ;a 1 p+1) p+1+O( p+1): Decomposing the function f into series
1 = ;(b 1 p+1 ;a 1 p+1) p ;(b 0 p+1 ;a 0 p+1)h p+1 1 +O(h p+2 1 ) +O( p+2) and neglecting the terms of order O(hp+1
1 ) and O( p+1), we obtain thatx(ti+ ) satises (46) dxd =f(x) + (p+ 1)(b1 p+1 ;a 1 p+1) p:
It can be shown by analogy that function x(ti+ ) on the every interval titi+1], i= 2:::N ;1 satises (47) dxd =f(x) + (p+ 1)(bi p+1 ;a i p+1) p where =(ti+1 ;t i), 2 01].
Thus the functionx(t) on 0T] (excluding grid nodes) satises (48) dxdt =f(x) +(t)
where function(t) has a piecewise-polynomial form (49) (t) = (p+ 1)(bi p+1 ;a i p+1)(t ;t i;1) p fort2 t
i;1ti] and i= 1:::N. Since on every subinterval ti;1ti] of the grid DN the equation (48) is satised with order of accuracy of O(hp+1
i ), we can replace (t) with a piecewise-constant function
'(t), which fort2 t i;1ti] equals (50) '(t) = (p+ 1) Z 1 0 (bi;1 p+1 ;a i;1 p+1) (t ;t i;1) pp d= = (bi;1 p+1 ;a i;1 p+1)h p i i= 1:::N: Consequently, instead of (48) we can consider
(51) dxdt =f(x) +'(t):
Recall that the sequence fD N
g of grids on 0T] was chosen in such a way that there exists a piecewise-constant function
: 0T]! C
1C2] 0< C1 C2 such that for every step of a gridDN, the following holds:
(52) hi =
T
N (ti;1) i= 1:::N
where N is the number of subintervals in the DN grid. Then the equation (51) can be written as
(53) dxdt =f(x) + (bp+1 ;a p+1)(t) T N (t) p where (bp+1 ;a p+1)(t) = (b i;1 p+1 ;a i;1 p+1) fort2 t i;1ti] and i= 1:::N.
Let us choose equal to (T=N)p. Then for bigN the global dis-cretization error can be asymptotically decomposed in a power series expansion (3) with respect to the small parameter and the main term is the solution of the initial value problem
(54) _ =f0 x(x(t))+ (b p+1 ;a p+1)(t) p(t) (0) = 0
being similar to (4).
Since the local error for method (24) is (55) x n+1 ;x(t n+1) = xn+1 ;x n+1+O(h p+2 n+1) = = (bn+1 ;a n+1)(tn) h p+1 n+1+O(h p+2 n+1) and p= 3, we can take
(56) Np Tph n+1 ( xn+1 ;x n+1) as the term (bn+1 ;a n+1)(tn) p(t
n) in the discrete problem.
4. Algorithm
To solve the initial value problem (1), we chose the initial gridDN 0 =t0 < t1 < ::: < tN =T hi=ti
;t
i;1 i= 1:::N and used the third-order nite-dierence scheme (see (24) or (28)):
(57) x n+1 ;x n hn+1 = 12 K1 ; 3 2 K3+ 2K4:
Moreover, at every step we calculated xn+1 using the fourth order R-K-M scheme (see (25) or (29)): (58) xn+1 ;x n hn+1 = 16 ( K1+ 4K4+K5) wheren= 01:::N;1 and K1 =f(xn) K2=f(xn+hn+1=3K1) K3 =f(xn+hn+1(K1+K2)=6) K4=f(xn+hn+1(K1+ 3K3)=8) K5=f(xn+hn+1 (1=2K1 ;3=2 K 3+ 2K4)):
We used (T=N)3 as a small parameterin system (2), where N is the number of subintervals inDN and (tx) is calculated as (59) (tn+1xn+1) = N3 T3h n+1 ( xn+1 ;x n+1) and xn+1 ;x n+1 is dened by (30) forn= 01:::N ;1.
To calculate the numerical solution of the initial value problem (4), we used modied Euler scheme with recalculation
n+1=n+hn+1 f0 x^+(t n+1xn+1) where f0 x= f 0 x(x n+1) +f 0 x(x n)]=2 ^ =n+hn+1=2 f 0 x n+(tn+1xn+1)] 0= 0:
5. Numerical experiment
We conducted numerical experiment for system 5,6] (60) 8 < : _ x1=k1(1 ;x 1 ;x 2) 2 ;k ;1x 2 1 ;2k 3(x2) x 2 1x 2 f 1 _ x2=k2(1 ;x 1 ;x 2) 2 ;k 4(x2x3) x 2 ;k 3(x2) x 2 1x 2 f 2 _ x3=" x2(1 ;x 3) ;x 3(1 ;x 1 ;x 2)] g where k3(x2) =k30exp( ; 3x2) k4(x2x3) =k40exp( ; 4x2 ; 5x3) with the parameters
k1= 0:2 k;1= 0:0025 k2 = 15 k30= 100 k40= 2
3 = 30 4 = 12 5 =
;10 = 40 "= 0:000025: The dynamic behavior of the system (60) is shown in the Fig. 1 for 0< t <800 with the initial data x1 = 0:385672,x2 = 0:5061437,
x3 = 0:4345207. The projections of the corresponding approximate solution on dierent phase planes (xixj) are shown in Fig. 1 as well. Using the described above algorithm we found the vector function
to be used in (2) for the R-K scheme (24) and the solution of the initial value problem (4), namely the vector function . Fig. 2 and 3 present the vector-functions,andand their projections on the dierent coordinate planes.
Note that our calculation reaches asymptotic behavior, i.e., that functiondoes not depend on when 10;8< <10;6. However, if we decreasedfurther (grid rening) the absolute value of the func-tion started increasing. The following explains why we did observe that previously 6].
To determine the function (59) we have to take a product of (61) x n+1 ;x n+1 hn+1 = 13 K1 ; 9 2K3+ 4K4 ; 1 2K5 :
and (N=T)3 = ;1, that is why during the calculation we have to provide following condition: (61) must be bigger then 10;11. In our case when it was < 10;9, the roundo errors became bigger then reminder term of the formula (57), i.e. (61) became smaller then 10;11and was computed wrong. That was proved by computations:
graph image looked like a "saw tooth curve".
However if the calculations were conducted with double precision, the eect disappeared and function did not change with decreasing of . Thus with the subtraction of the quantities of the same order
of magnitude in (60), we had seen a loss of signicant digits that did lead to errors in computation of.
To decrease the addition and subtraction error in (60), we present the system (60) as _
x
=P
(x
) +Q
(x
) wherex
= (x1x2x3) TP
= (f 1f2g) T ;Q
Q
= (;2k 3(x2)x 2 1x 2 ;k 3(x2)x 2 1x 2 0) T: Because on the solutions that we are interested in there isk
Q
(x
)k kP
(x
)kand vector
Q
strongly changes with time, we calculated the right hand side of (60) separately for vectorsP
andQ
, and after that added them together.References
1. Dekker, K. and Verwer, J.G. (1984): Stability of Runge-Kutta methods for sti nonlinear dierential equations,North-Holland, Amasterdam{New York{ Oxford.
2. Godunov, S.K. and Ryaben'kii, V.S. (1973):Finite-dierenceschemes, Nauka, Moscow. English Transl.:Dierence schemes. An introduction to the underly-ing theory, Studies in Mathematics and its Applications, 19, North-Holland,
Amsterdam (1987).
3. Stetter H.J. (1973):Analysis of Discretization Methods for Ordinary Dier-ential Equations, Springer-Verlag, New York.
4. Chumakov, G.A., Slinko, M.M., Belyaev, V.D., and Slinko, M.G. (1977): Ki-netic model of an oscillating heterogeneous reaction, Dokl. Akad. Nauk SSSR
234, No. 2, 399{402.
5. Chumakov, G.A. (1980): Mathematical modeling of complexe behaviour of heterogeneous catalytic reaction rate, Report No. 233, Computing Center SB AS USSR, Novosibirsk.
6. Chumakov, G.A. and Slinko, M.G. (1982): Kinetic turbulence (chaos) of re-actions rate for hydrogen oxidation on metallic catalysts, Dokl. Akad. Nauk SSSR266, No. 5, 1194{1198.
7. Kubicek, M. and Marek, M. (1983):Computational Methods in Bifurcation Theory and Dissipative Structures, Springer-Verlag, New York.
8. Guckenheimer, J. and Holmes, Ph. (1997):Nonlinear Oscillations, Dynamical Systems and Bifurcation of Vector Fields, 5th Edition, Springer-Verlag, New York.
9. Cano, B. and Sanz-Serna, J.M. (1997): Error growth in the numerical in-tegration of periodic orbits, with application to Hamiltonian and reversible systems, SIAM J. Numer. Anal.34, No. 4, 1391{1417.
Fig.1. Dynamics of the kinetic model (60) and their phase plots. Here X, Y and Z correspond tox 1, x 2, and x 3, respectively
Fig. 2. Dynamics of the local and global errors of the kinetic model (60). Here X, Y and Z correspond tox 1, x 2and x 3, respectively
Fig.3. Local and global errors in dierent phase-planes for the system (60). Here X, Y and Z correspond tox 1, x 2and x 3, respectively
