CC BY
30
Fig. 3. The noise immunity of IRAM at k = 0,7 and INV, = 1; INV. = 2; 3; 4; 5; 6; 7.
l 5 i 5 5 5 5 5
Curve 1is the probability of the pairwise transition into IRAM under non-correlativity of the noise readings (theoretical limit; Curve 2 is the probability of the pairwise transition of IRAM at R = 0,7; Curve 3 is the error probability of the classical AM
Bibliography
1. The Invariant Method of Analog ofTelecommunication System of Information Transmission : monograph /
V. B. Malinkin, E. I. Algazin, A. N. Levin, N. V. Popantonopulo. Krasnoyarsk, 2006.
2. Algazin, E. I. The Estimation of the Noise Immunity of the Invariant System ofInformationProcessing underNon-Coherent Reception / E. I. Algazin, A. P. Kovalevsky,
V B. Malinkin. Vestnik SibSAU, Krasnoyarsk, 2008. Iss. 2(19). P. 38-41.
3. Algazin, E. I. The Noise Immunity of the InvariantRelative Amplitude Modulation / E. I. Algazin, A. P. Kovalevsky,
V B. Malinkin//MaterialsofIXIntem. Conf. “ActualProblems of Electronic Instrument Engineering” (APEIE-2008). Novosibirsk, 2008. Vol. 4. P. 20-24.
4. Algazin, E. I. The Noise Immunity of the Invariant System of Information Transmission under the Weak Correlation Communications /E.I. Algazin, A. P. Kovalevsky, V. B. Malinkin // Vestnik of SibSAU. Krasnoyarsk,2008, P. 29-32.
5. Teplov, N. L. The Noise Immunity of the Systems of TransmissionofDiscrete Information/N.LTeplov. M., 1964.
6. Borovkov, A. A. The Theory of Probabilities / A. A. Borovkov. M., 1999.
7. Levin, B. P. Theoretical Foundations of Statistic Radio Engineering/B. P. Levin. M., 1989.
© Algazin, E. I., Kasatkina E. G., KovalevskyA. P., Malinkin V. B., 2009
A. E. Novikov SiberianFederal University, Russia, Krasnoyarsk
E. A. Novikov
Institute of Computational Modelling, RussianAcademy of Sciences, Siberian Branch, Russia, Krasnoyarsk
THE NUMERICAL MODELING OF A CESIUM CYCLE IN THE UPPER ATMOSPHERE BYAN L-STABLE METHOD OF SECOND-ORDERACCURACY*
An algorithm of right-hand side and Jacobian formation of differential equations of chemical kinetics is described. Numerical simulation of the cesium cycle in the upper atmosphere is conducted by means of the L-stable method of the second order of accuracy with the control accuracy. The results of the computation are presented.
Keywords: chemical kinetics, cesium cycle, stiff problem, L-stable method, accuracy control.
The modeling of chemical reactions kinetics’ is applied in the studies of various chemical processes. The subject of this study is the time dependence on concentration of reagents being a solution for the Cauchy problem and for systems of ordinary differential equations. Difficulties in solving such problems are related to stiffness and a large scale. In modern methods forthe solving stiff problems, an inversion of the Jacobi matrixoftheequations’ systemisused. Inthecaseoftheoriginal problem’s large scale, the decomposition of the given matrix essentially defines a total of computational work. To improve calculation efficiency, in a number of algorithms the freezing of the Jacobi matrix, i.e.the application of the same matrix at several
iteration steps, is used [1-2]. This approach is used in advantage to algorithms based on the multistep methods, in particular the formulaforbackwarddifferentiation [3]. The situationis worse in algorithms for integrations based on known iteration-free methods among which are the Rozenbrok type methods and theirmodifications [1].Hers isanalgorithmforthe construction of the right-hand side and the Jacobi matrix of the differential chemical kinetics’ equations. Results of numerical modeling of an ionization-deionization cesium cycle in the upper atmosphere withan L-stable method of second-order accuracy, in which the numerical freezing as well as analytical Jacobi matrix is allowed, aregivenhere.
* The work was supported by the Russian Foundation of Fundamental Researches (grant №08-01-00621) and by the Presidential grantNSh-3431.2008.9.
A kinetic scheme of a chemical reaction consisting of elementary stages which are written in formula [4]
Nr kj Nr
(aijxi----------------------------------k(P ijX1 (1)
i=1 i = 1
wherex, 1<I<Nr are reagents; k, 1<j<Ns are constants of stages’ rates; Nr and Ns are a number of reagents and a number of stages respectively; a, and Mif 1 < I < Nr, 1 < j < Ns are stoichiometric coefficients. Inthe framework of an isothermal reactor of a constant volume lumped model, to the process (1) corresponds a system of ordinary differential equations C ’ = ATV (2)
with a given initial condition C(0) = C0. Here AT is a stoichiometric matrix, C and V are vectors of reagent concentrations and rates of stages respectively. When a reaction occurs in non-isothermal conditions, the heat balance equation is
T ' = [QTV -a (T - T,,)] / CTVC,
(3)
where T is the temperature of a reactor mixture, T01 is temperature of reactor walls, QT is a vector of specific heat of stages, CTV is avector of the reagents’ heat capacity, a = a s/ r, a is a heat transfer coefficient, s and r are an area of a surface and the volume of a reactor respectively, is added to the system (2). The superscript T of the vectors QT and CTV means transposition. Heat capacity of reagents and a heat transfer coefficient may be the functions of reagent concentration a may depend on temperature as well.
When the reaction occurs in an isothermal reactor of a constant volume with an agent change (an open system, a reactor of the ideal mixing), a system of differential equations is written in the formula
C' = AVT +±(Cp -C),
(4)
K= ks n c“" i W-= k- s n cf“, i=1 i=1
where k and k ,1 < s < N are constants of rates of direct and
s -s7 s
reverse stages, respectively. Constants of rates are calculated by the formula
kj = AjTnj exp(-Ej / RT), where T is temperature of the mixture in the reactor; A, nj and EJR are the given constants. The immediate use of this formula may lead to a wrong result or arithmetic overflow because of the large constant data spread [5-6]. Therefore for calculating constants of stages’ rates the following relation is used:
kj = exp(ln Aj % nj ln T - Ej / RT).
The stoichiometric matrix AT with entries a., is formed from the kinetic scheme (1) by the following rule: a stage number coincides with a column number and a reagent number coincides with a number of the row of the matrix AT. If x. is an initial reagent then a = -ajf if x. is a product then a = M If x. is an initial reagent and a product simultaneously then a.. = My — a... Usually a small number of reagents take part in an elementary stage, i. e. a stoichiometric matrix is sparse. For the system (2) the j-th column b of the Jacobi matrix defined by the formula
T RV
bj = AT — i 1 < j < Nr.
oc,
(6)
where Cp is a vector of reagent concentration in the inlet, © is time of the mixture being in the reactor, © = r/u, u is volume velocity of the mixture flow through a reactor. When a reaction occurs in non-isothermal conditions, the heat balance equation is
T' = [QTV -a (T - Tm)] / CTVC -1(T - ^), (5)
where T02 is the temperature of a mixture at the inlet of a reactor, is added to the system (4). The temperature of a reacting mixture can be given as a function of time and concentration, i. e. T = T(t,C).
If a stage is reversible, then the difference of rates W+ and W— (of direct and inverse processes) respectively, is conventionally taken as a rate Ws of a stage, i. e.
ws = w%% w-, 1 < s < ns .
If the third particle is involved in the reaction, then a rate is recalculated by the formula
Nr % Nin
Vs = PsWs 1 Ps = ( Ssict, i =1
where e , 1 < i < N are efficiencies of the third particles, N ,
sv r L 7 in7
e and c, N +1 < i < N + N are quantity, efficiencies and
si v r r in A J 7
concentrations of inert agents, respectively. The values of the components of the vector Ws are determined from scheme (1)ofa chemical reactionby the relations
In the case that the system (2), (3) is solved simultaneously, the row bNr + ip the column bN +j and the corner element bNr + + p which are determined the following
way:
bK = !QT RV(T - t„) %
r oc. oc.
% [a(T - T») - QTV ]cj CTC} / CTC,
bi ,Nr %1 = (A -RT )i1 1 < i < Nr 1
r oT
t RV Ra t
bNr %1 Nr+1 = [QT—(T - Tn) ^ ] / CV -
r r RT RT
- [QT V -a (T - T)1)](C) / (CT C )2,
RT
are added to the matrix. For a flow reactor, the diagonal elements of the obtained matrix are adjusted considering the component 1/©. Observe that the components of the vector
RV/Rc have the form
i
Nr % Nin
c^s /Rc< =«siPskscT'~l n c<st -
-PsiPsk-scfs-1 n cisk -^(W% - w-),
(7)
and for the components of the vector RV/RT we have the following reactions:
Rk Nr % Nin RVs /RT = Ps n c“” -
Rk
RT
n c?" i 1 < s < Ns.
If at the s-th stage there is no third particle, then in the expressions for Rv /Rc and Rv/RT one should put p =1 and e =0.
s i s s si
To determine the entries b , 1 < i, j < N of the Jacobi
if 7 J r
matrix withthe help of (6), we can apply the formula
b> as ft' 1 < '•j < N'-
Considerthe individualterm(7), i. e.
Nr & Ntn
a. p k a .ca"-1 ) ca'k -
jsr^s s si i X X i
k =1#k (i
Nr & Nn
-a . pk B cf”-1 ) c^’k & as. (W + - W~).
jsr s -sr si i X X i js si \ s s J
k=1# k (i
To determine this expression, it is necessary to perform three steps. Repeating them in a loop allows the forming of the relations (8). At the first step apk a and apk + , are
v ' A js s s si js s —s si7
formed with pk or pk assumed to be calculated. At the
r s s r s —s
second step we define
ae
c““-1 n c°'k,
1 n
In the case that the third particles are involved in the scheme of a reaction, for the third step a es( W+- W-) is given.
We describe a numerical method which is applied for numerical integration of an ionization-deionization cesium cycle inthe upperatmosphere. Forthe numerical solution of the Cauchy problem in a system of ordinary differential equations
y' = f (t,y), y(to) = yo, t0 < t < tk, (8)
where y and f are real N-dimensional vector functions, t isan independent variable, we consider the (m, k)-method of the form [7]
V2
y„&1 = yn & ak1 &(1 - a)k2' a =1 -—,
Dnk1 = hnf (tn &Ph' y„)' Dnk2 = k1. (9)
Here k. and k. are stages of the method; D = E- ah A ;
12° 7 n n n7
E is the identity matrix; hn is an integration step; Anisamatrix whichcanberepresentedintheform A = f +hB + O(h2);
n n n
f n = ff(yn)/dy is the Jacobi matrix of the system (8); Bn is a matrix which does not depend on an integration step. We can apply the method (9) of freezing the numerical as well as the analytical Jacobi matrix. Using the Taylor expansion of the stages, alocalerror 8n ofthemethod(9)canbewritteninformula
S„ = (a -1 /3)h3fy7&h3[^4 f-&6 &
& 3 ftf - 22 ff - ^2 Bf ] & 0(h4).
Then, according to [8], forthe accuracy control in (9) we
can use the error estimation of the form
e„ = (a -3)h2f;f & O(h3).
Taking into account
k2 - k1 = ah f'y,„f„ & O(h3)
we can estimate en within the order terms of O(h3) using this formula
e:„ = a-1(a -3)[k2 - k^.
As a result, in the accuracy control for the method (9)we can apply the inequality
e(j„) =11 D„" k -*,)||< 1/3 .
1/3-a
where e is a required accuracy of integration, ||-|| is some norm in RN, and the value of the parameter j, 1 < j < 2, is
7 i J n7 J n 7
taken the least one for which this inequality holds. In precise calculations the norm in the inequality had been calculated withtheformula
IID1-j„ k )ll |[D-jn (k2 -k,)]i |
11 Dn (k2-k')||= !sax—iyr&r— ’
where i is the number of the component of the approximate solution, tr is a positive parameter. If for the i-th component of a solution the inequality \y‘n| < tr holds then the absolute error tr-e is controlled, inthe opposite case the relative error e is controlled. Inthe calculations the parameter tr is chosen so that for all components of the solution the actual accuracy is not worse than the given one.
In the numerical calculation of the Jacobi matrix a step r 1 < j < N of numerical differentiation is takenby the formula r. = max(10~14,107yj|). Usually inthe solution of stiff problems double precision is used since the Jacobi matrix of the system (9) is ill-conditioned. Now the j-th column aj of the matrix A
nn
in (9) is calculated by the formula
f(ylyj & rjyN) - f(ylyjyN)
an = ----------------------------------------'
rj
i. e., calculation of the matrix An requires N calculations of the right-hand side of the problem (8). An attempt to use the previous matrix Dn is performed after each successful integration step. To ensure that the stability property of a numerical scheme holds, the integration step remains constant when freezing the matrix Dn. U nfreezing of the matrix takes place in the following cases:
- the accuracy of calculations is violated;
- the number of steps with a frozen matrix achieves a given maximal value Ih,
- the predicted step is Qh times greater than the last successful one;
- e(1)>e(2).
The choice of Ih and Qh can have influence on redistribution of computational work. For Ih =0 and Qh = 0 freezing a matrix does not take place, with increasingIh and Qh the number of calculations of the right-hand side of the problem (8) increases and the number of inversions of the Jacobi matrix decreases. Efficiency of the integration algorithm depends on Ih and Qh by the choice of which the algorithm can be adjusted for a concrete type of problems. They are chosen depending on the ratio of the computational cost of the function f to that of the decomposition of the Jacobi matrix.
Now we consider a model of the ionization-deionization cesium cycle in the upper atmosphere. The presented model is obtained from a large kinetic scheme and is widely used abroad as a typical example of a stiff problem ofkinetics [9]. The scheme of the reaction has the form
1) O&&C& 1 Cs &O2,
2) C& & e -> Cs,
3) Cs 1 C& & e,
4) O2 & Cs & M 1 CsO2 & M,
5) O2 & e & M 1 O2 & M ,
6) O2 102 & e,
where the constants of rates are the following: k1 = 3,0-1010,
k2= 6,0-105, k3= 3,24-10-3, k4 = 3,63-104, k5 = 3,63-104, k6= 4,0-10-1. The reaction occurs with participation of the inert substance N2, its concentration is [N2] = 3,32-10—3. Efficiencies of the third substances are equal to 1 for all reagents with the exception of the efficiency of O2 at the fifth stage which is equal to 12,4. We denote the concentrations of the reagents in the following way:
c1 = [eL c2 = [O2 L c3 = [Cs ],
c4 = [CsO2], c5 = [Cs+ ], c6 = [O2].
The corresponding system involves 6 differential equations of the form
c1 = k2c!c5 % k3c3 - k5P2c1c6 % k6c2 ,
c2 = k1 c2 c5 % k5 p2 c1c6 - k6 c2 ,
c3 = klc2c5 % k2clc5 - k3c3 - k4pc3c6,
c4 = k4 pxc3c6, (10)
c5 = —k^c5 — k2cc % k3c3,
c6 = k1c2c5 - k4p1c3c6 - k5p2c1c6 % k6c2,
where
p1 = c1 % c2 % c3 % c4 % c5 % c6 % [N2 ] ,
p2 = c %c2 %c3 %c4 %c5 % 12.4c6 %[N2].
Integrationisperformedoverthesegment [0.1000] with initial step 10-5 for the following initial concentrations of the reagents:
c1 = [e] = 1.66 -10 16, c2 = [O2- ] = 8.63 -10 16, c3 = [Cs] = 1.66-10 6, c4 = [CsO2] = 0.0, c5 = [Cs+ ] = 1.03 -10-15, c6 = [O2] = 5.98-10 4.
The presented algorithm was compared in efficiency with the known dlsode method of Gear in the implementation of A. Hindmarsh [10] foraccuracy e = 10-2 of calculations. Poor accuracy of calculations is related to the fact that the method (9) is of second-order accuracy and there is no sense in using it in high-accuracy calculations. The number if of calculations of the right-hand side and the number ij of decompositions of the Jacobi matrix of the problem (10) on the integration interval are taken as an efficiency criterion. To solve the problem (10), the constructed algorithm requires
101 calculations of the right-hand side and 14 decompositions of the Jacobi matrix. In the dlsode method the required accuracy of 10-2 is achieved for e = 10—3 with computational costs if= 192 and ij = 22.
The numerical results show that the proposed algorithm has the advantage for the number of calculations of the right-hand side almost by a factor of two, whereas the advantage for the number of decompositions of the Jacobi matrix is aboutby a factor of1.5. This integration algorithm is assumed to have the most efficient applicationfor stiff problems for accuracy e=10-2 of calculations (of order of one percent) and lower.
Bibliography
1. Hairer, E. Solving Ordinary Differential Equations. Stiff and Differential — Algebraic Problems / E. Hairer, G. Wanner. Springer-Verlag, 1996.
2. Novikov, E. A. Explicit methods for stiff systems / E. A. Novikov. Novosibirsk: Nauka, 1997.197 p. (inRussian).
3. Byrne, G. D. ODE solvers: a review of current and coming attractions / G. D. Byrne, A. C. Hindmarsh // J. of Comput.Physics. 1987. № 70.P 1-62.
4. Emmanuel, N. M. Chemical Kinetics / N. M. Emmanuel,
D. G. Knorre. Moscow : Vysshaya shkola, 1974. 400 p. (in Russian).
5. Babusok, V. I. Numerical solution of direct kinetic problems/V. I. Babusok, E. A. Novikov //React. Kinet. Catal. Lett. 1982.V21.№. 1-2.P. 121-124.
6. Babusok, V. I. Generator of a right-hand side and a Jacobi matrix of differential equations of chemical kinetics : preprint № 359/V. I. Babusok, E. A. Novikov ; Computing Centerof SB AS USSR. Novosibirsk, 1982.27 p. (inRussian).
7. Novikov, E. A. The (2,1)-method for stiff non-autonomous problems/E. A. Novikov // Systems of control and informational technologies. 2008. № 2(32). P. 12-15. (in Russian).
8. Novikov, E. A. Some methods for stiff systems induced by one and two calculations of a right-hand side /
E. A. Novikov, Yu. A. Shitov // Math. models and methods forproblems of mechanics of continua. Krasnoyarsk, 1986. P. 11-19.
9. Edelson, D. The new book in chemical kinetics / D. Edelson//J. Chem. Ed. 1975.V 52. P. 642-643.
10. Brown, P. N. Reduced Storage Matrix Methods in Stiff ODE Systems / P. N. Brown, A. C. Hindmarsh // J. Appl. Math. & Comp. 1989. V31.P. 40-91.
© Novikov, A. E., Novikov, E. A., 2009
