- Email: [email protected]
Explicit appropriate basis function method for numerical solution of stiff systems
Annals of Nuclear Energy 75 (2015) 353–357
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Explicit appropriate basis function method for numerical solution of stiff systems Wenzhen Chen a,b,⇑, Hongguang Xiao a, Haofeng Li a, Ling Chen a a b
Department of Nuclear Energy Science and Engineering, Naval University of Engineering, Wuhan 430033, China Institute of Thermal Science and Power Engineering, Naval University of Engineering, Wuhan 430033, China
a r t i c l e
i n f o
Article history: Received 11 April 2014 Received in revised form 7 July 2014 Accepted 20 August 2014
Keywords: Stiff systems Basis function Ordinary differential equations Explicit method
a b s t r a c t In this paper, an explicit numerical method, called the appropriate basis function method, is presented. The explicit appropriate basis function method differs from the power series method because it employs an appropriate basis function such as the exponential function, or periodic function, other than a polynomial, to obtain approximate numerical solutions. The method is successful and effective for the numerical solution of the first order ordinary differential equations. Two examples are presented to show the ability of the method for dealing with linear and nonlinear systems of differential equations. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction Numerous works have been focusing on the development of more advanced and efficient methods for solving ordinary differential equations (ODEs) (Butcher, 2000; Hajmohammadi et al., 2012, 2014a). For example, Wazwaz (2010) presented a new algorithm for solving in ODE’s of integro–differential type, and Hajmohammadi and Nourazar (2014b, c) presented a new algorithm for solving in ODE’s of eigen-value type based on semi-analytical method and pure-analytical method for solving in ODE’s of non-linear type, respectively. Implicit numerical methods are usually used for solving ordinary differential equations, in particular for some stiff problems, in which the implicitness will cause the method very involved compared to explicit methods (Li et al., 2009; Zhu et al., 2012). Indeed, explicit numerical methods require smaller step sizes in situations where implicit methods would not. It should be emphasized that explicit methods also show important advantages, some authors have developed many explicit methods for stiff systems (Guzel and Bayram, 2005). There has been a great deal of research that focused on eliminating the stiffness problem of reactor kinetics (Koclas et al., 1996; Chen et al., 2013). And there are several methods especially adapted for solving the initial value problems for stiff systems of ⇑ Corresponding author at: Department of Nuclear Energy Science and Engineering, Naval University of Engineering, Wuhan 430033, China. Tel.: +86 027 13871167436. E-mail address: [email protected] (W.Z. Chen). http://dx.doi.org/10.1016/j.anucene.2014.08.040 0306-4549/Ó 2014 Elsevier Ltd. All rights reserved.
ordinary differential equations (Aboanber and Hamada, 2003; Aboanber, 2004; Tashakor et al., 2010). Among the methods are numerical integration using Simpson’s rule, finite element method, Runge–Kutta procedures, quasi-static method, piecewise polynomial approach and other methods (Li et al., 2010; Abdallah and Nahla, 2011; Hamada, 2013). Most of these methods are successful in some specific problems, but still suffer, more or less, from disadvantages as mentioned by Chen et al., 2013. In this paper, an explicit numerical method for stiff systems is developed and tested. The method differs from the power series method (Guzel and Bayram, 2005) in the idea that employs an appropriate basis function, such as the exponential function, periodic function, and so on, other than a polynomial to obtain approximate numerical solutions. In the other hand, the method is considerably more accurate than the power series method, as demonstrated in following sections. 2. Appropriate basis function method 2.1. Basic ideas Since every ordinary differential equation with order n can be written as a system consisting of n ordinary differential equations with order one, our study is restricted to a system of first order differential equations, which can be written as follows
y0 ¼ f ðt; yÞ; t 2 ½0; T; kyk < þ1; yð0Þ ¼ y0 ;
ð1Þ
354
W. Chen et al. / Annals of Nuclear Energy 75 (2015) 353–357
whose theoretical solution is y(t). Let yn be an approximation to y(tn). Here
terms, Eq. (4) can also give the very accurate approximation. For example, when m = 2, Eq. (4) can be written as follows
y ¼ ½y1 ; y2 ; . . . ; yk T ;
tanhðtÞ 1 2e2t þ 2e4t e5t ;
T
f ¼ ½f 1 ; f 2 ; . . . ; f k ; and
y0 ¼ ½y01 ; y02 ; . . . ; y0k T : It is assumed that f and y are sufficient differentiability (Guzel and Bayram, 2005). The power series method (Guzel and Bayram, 2005) assumed that the solution of Eq. (1) can be expressed by a polynomial expression as follows
y¼
1 X Ai xi ;
ð2Þ
i¼0
where Ai is a vector function which is the same size as y0. However, the authors think that it should choose an appropriate function as the basis function according to the characteristic of stiff systems in order to improve the computational efficiency and accuracy. It is well known that a real function can be expressed by some different basis functions. For example, let us consider u(t) = tanh (t), which can be expressed by a polynomial expression as follows
tanhðtÞ t
m!þ1
1 X Ai xi ðxÞ;
ð5Þ
i¼0
where xi(x) is the ith term of the basis functions. For example, if the set of appropriate basis functions is fexpðixÞji 0g , xi(x) will be exp (ix).
ð3Þ
Then Eq. (3) converges to the exact solution u(t) only in a small region 0 6 t < 0.5 (see Fig. 1). However, by means of the exponential function as the basis function, one has
tanhðtÞ 1 þ lim 2
y¼
2.2. Scheme of the method
t 3 2t 5 17t 7 þ þ ; 3 15 315
"
which agrees well with the exact solution tanh (t) in the whole region 0 6 t < + 1. In addition, the numerical tests in Table 1 demonstrate that the efficiency of the exponential function is better than that of a polynomial expression. When the approximation solution is expressed by the exponential function u(t), it will take a computing time of 0.094 s to obtain the value of u(t) in the time interval [0, 3/2] with the time step h = 0.00001 s. However, using a polynomial as the basis function, it will take a computing time of 0.203 to 0.297 s. Therefore, one can get the best approximation by means of an appropriate basis function. For example, a periodic solution is expressed more efficiently by periodic basis functions than by a polynomial expression. The solution of Eq. (1) can be expressed by a set of appropriate basis functions fxi ðiÞji 0g as follows
#
m X
ð1Þn e2nt þ ð1Þmþ1 eð2mþ1Þt ;
The way that the appropriate basis function method is used in practice is carried out by the following computation steps: Step 1: Choosing the order of approximation solution, m. Step 2: When x = 0, according to Eq. (5), we can get
ð4Þ
yðnÞ ð0Þ ¼
n¼1
which converges to the exact solution u(t) in the whole region 0 6 t < + 1 (Liao and Tan, 2007). And, even taking the first few
m X ðnÞ Ai xi ð0Þ:
ð6Þ
i¼0
Step 3: According to Eq. (1) and the initial conditions, the values or expressions of yð0Þ; y0 ð0Þ; y00 ð0Þ; . . . ; yðmÞ ð0Þ can be obtained by the iterative computation. Step 4: When n ¼ 0; 1; 2; . . . ; m, according to Eq. (6), the linear equations of Ai can be derived, and it is easy for us to develop a computer code to calculate the coefficient Ai. Then, substituting Ai into Eq. (5), the m th-order approximation solution of Eq. (1) is obtained . Step 5: Making step size of x to be h and substituting it into the mth-order approximation solution of Eq. (1), we have y at x = x0 + h. Step 6: Repeating above steps 1 to 5, we can obtain the numerical solution of Eq. (1). 3. Numerical experiments The appropriate basis function method and the power series method (Guzel and Bayram, 2005) are compared in the following cases. We arranged the independent variable x, the step size h,
Fig. 1. Exact and approximation solutions of u(t) = tanh (t).
Table 1 Comparison of approximation solutions of u(t) = tanh (t) in the time interval [0, 3/2] by using h = 105. Type of approximation
Solutions for different time t = 0.8
t=1
Computing time (s) t = 1.1
t = 1.2
t = 1.3
t = 1.4
t = 1.5
uðtÞ ¼ t t3 þ 2t5 17t 315
3
5
0.7491
1.0127
1.1954
1.4260
1.7142
2.0677
2.4904
0.297
t3 3
5 þ 2t5 2t
0.7604
1.0667
1.3005
1.6193
2.0528
2.6366
3.4125
0.203
0.6594 t = 0.8 0.6640
0.7592 t=1 0.7616
0.7989 t = 1.1 0.8005
0.8326 t = 1.2 0.8337
0.8610 t = 1.3 0.8617
0.8849 t = 1.4 0.8854
0.9048 t = 1.5 0.9052
0.094
uðtÞ ¼ t
7
þ 2e4t e5t uðtÞ ¼ 1 2e Accurate solution u(t) = tanh (t)
355
W. Chen et al. / Annals of Nuclear Energy 75 (2015) 353–357 Table 2 Comparison of theoretical and numerical results in problem 1. Numerical results h
N
y
Theoretical results
Appropriate basis function method (i = 4)
1 ¼ 1 0.05
1
1
1
y1 y2 y1 y2
0.90483741803596 0.95122942450071 0.13533528323661 0.36787944117144
Power series method (Guzel and Bayram, 2005) (i = 4)
1 ¼ 0:6
Solution
Error
Solution
Error
Solution
Error
0.90483741803596 0.95122942450071 0.13533528323661 0.36787944117144
0 0 0 0
0.90483743957714 0.95122944625838 0.13611513643068 0.36883051686970
2.1541e8 2.1758e8 7.7985e4 9.5108e4
0.90483750000000 0.95122942708333 0.33333333333333 0.37500000000000
8.1964e8 2.5826e9 0.197998 0.007121
the computation step number N, the theoretical solution, and the numerical results by the appropriate basis function method and the power series method, respectively. All computations are carried out by the MATLAB code with double precision arithmetic.
y001 ð0Þ ¼ 1002y01 ð0Þ þ 2000y2 ð0Þy02 ð0Þ ¼ 4 y002 ð0Þ ¼ y01 ð0Þ ð1 þ 2y2 ð0ÞÞy02 ð0Þ ¼ 1
00 02 00 y000 1 ð0Þ ¼ 1002y1 ð0Þ þ 2000y2 ð0Þ þ 2000y2 ð0Þy2 ð0Þ ¼ 8 00 00 02 y000 2 ð0Þ ¼ y1 ð0Þ ð1 þ 2y2 ð0ÞÞy2 ð0Þ 2y2 ð0Þ ¼ 1
3.1. Problem 1
ð16Þ
;
:
ð17Þ
Then, the linear equation can be given as follows The test problem considers the following system of differential equations
where
y01 ¼ 1002y1 þ 1000y22
with initial conditions: y1(0) = 1 and y2(0) = 1. The exact solutions are y1(x) = e2x and y2(x) = ex. y1 and y2 can be represented by the exponential basis function other than a polynomial as follows, respectively 1 X
a1;i e1ix and y2 ðxÞ ¼
i¼1
1 X a2;i e1ix ;
ð8Þ
i¼1
where 1 is a constant, it can provide us with great freedom to choose appropriate basis functions to approximate the stiff system. Then, the n th derivatives of y1 and y2 are given by 1 1 X X n n ðnÞ ðnÞ y1 ðxÞ ¼ a1;i ð1iÞ e1ix and y2 ðxÞ ¼ a2;i ð1iÞ e1ix : i¼1
ð9Þ
i¼1
When m = 4, the 4th-order approximation solutions of Eq. (7) are given by
y1 ðxÞ ¼
4 X
a1;i e1ix
4 X and y2 ðxÞ ¼ a2;i e1ix ;
i¼1
2
ð7Þ
y02 ¼ y1 y2 ð1 þ y2 Þ;
y1 ðxÞ ¼
Aa ¼ B;
6 60 6 6 61 6 6 60 6 A¼6 6 12 6 6 60 6 6 6 13 4 0
i¼1
ðnÞ
i¼1
y2 ð0Þ ¼ 1
y02 ð0Þ
¼ y1 ð0Þ y2 ð0Þð1 þ y2 ð0ÞÞ ¼ 1
0
0
1
1
1
1
1
1
0
0
0
0
0
0
1
1
1
2
2
2
0
0
0
0
2
1
2
1
12
13 13 13
0
0
0
0
13 13 13 13
1
0
0
0
1
0
3
7 17 7 7 07 7 7 17 7 7; 07 7 7 27 1 7 7 07 5
19 2
1 2
ð18Þ
1 41 1 ; 4 þ 71 1
2
3
1
Finally, the 4th-order approximation solutions are obtained. Note that, when 1 ¼ 1, the exact solutions y1(x) = e2x and y2(x) = ex can also be obtained. The numerical results are illustrated in Table 2. 3.2. Problem 2
ð13Þ
Let us consider the following system of differential equations
i¼1
y01 ð0Þ ¼ 1002y1 ð0Þ þ 1000y22 ð0Þ ¼ 2
0
1
0
þ 72 12 þ 12 13 ; 1 11 11 12 16 13 T : 6
ð11Þ
ð12Þ
ð14Þ
;
0
1612 413 ;4 þ 1411 þ 1412 þ 413 ;
6
And according to Eq. (7) and the initial conditions, following equations can be derived
y1 ð0Þ ¼ 1
0
a ¼ 4 þ 26 11 þ 612 þ 43 13 ; 6 1911 3
i¼1
4 4 X X n n ðnÞ a1;i ð1iÞ and y2 ð0Þ ¼ a2;i ð1iÞ ;
0
11 412 43 13 ; 4 þ 133 11 þ 23 12 þ 16 13 ; 1 11 3
y1 ð0Þ ¼ a1;1 þ a1;2 þ a1;3 þ a1;4 and y2 ð0Þ
y1 ð0Þ ¼
1
Solving this linear equation results in
ð10Þ
Substituting x = 0 into Eqs. (10) and (11) results in
¼ a2;1 þ a2;2 þ a2;3 þ a2;4 ;
1
a ¼ ½a1;1 ; a1;2 ; a1;3 ; a1;4 ; a2;1 ; a2;2 ; a2;3 ; a2;4 T :
i¼1
4 4 X X n n ðnÞ ¼ a1;i ð1iÞ e1ix and y2 ðxÞ ¼ a2;i ð1iÞ e1ix ;
1
B ¼ ½1; 1; 2; 1; 4; 1; 8; 1T ;
According Eq. (10), following equations can be gotten ðnÞ y1 ðxÞ
1
y01 y02 y03
¼ 20y1 0:25y2 19:75y3 ; y1 ð0Þ ¼ 1; ¼ 20y1 20:25y2 þ 0:25y3 ; ¼ 20y1 19:75y2 0:25y3 ;
y2 ð0Þ ¼ 0; y3 ð0Þ ¼ 1:
ð19Þ
The analytical solutions of Eq. (19) are
y1 ¼ ½e1=2x þ e20x ðcosð20xÞ þ sinð20xÞÞ=2; ;
ð15Þ
y2 ¼ ½e1=2x e20x ðcosð20xÞ sinð20xÞÞ=2; y3 ¼ ½e1=2x þ e20x ðcosð20xÞ sinð20xÞÞ=2:
ð20Þ
356
W. Chen et al. / Annals of Nuclear Energy 75 (2015) 353–357
Table 3 Comparison of theoretical and numerical results in problem 2. Numerical results h
N
0.01
20
0.1
2
0.025
y
y1 y2 y3 y1 y2 y3 y1 y2 y3 y1 y2 y3
20
0.1
Theoretical results
5
0.43950209815007 0.45147399867173 0.45336341936422 0.43950209815007 0.45147399867173 0.45336341936422 0.38936899538135 0.38940708916983 0.38939369390157 0.38936899538135 0.38940708916983 0.38939369390157
Appropriate basis function method (m = 6)
Power series method (Guzel and Bayram, 2005) (i = 6)
Solution
Error
Solution
Error
0.43950209807394 0.45147411086509 0.45336330759892 0.45326849690173 0.44691922771744 0.45791817689314 0.38936955757924 0.38940855772086 0.38939222627783 0.38939761019105 0.38939766195830 0.38940310409674
7.61273e011 1.12193e007 1.11765e007 0.013766398 0.004554770 0.004554757 5.62197e007 1.46855e006 1.46762e006 2.86148e005 9.42721e006 9.41019e006
0.43950209127250 0.45147399307674 0.45336342495922 0.42698661025269 0.51933228926504 0.38550512877121 0.38936730990046 0.38940763611512 0.38939314695629 0.39162253788905 0.38878523056801 0.39001555250403
6.87757e009 5.59499e009 5.59499e009 0.012515487 0.067858290 0.067858290 1.68548e006 5.46945e007 5.46945e007 0.002253542 6.21858e004 6.21858e004
In addition, Eq. (19) can be rewritten by a matrix formulation as follows
dy ¼ Hy; dx
ð21Þ
yðxÞjx¼0 ¼ y0 ;
ð22Þ
where y = [y1, y2, y3]T and y0 = [1, 0, 1]T. H is a square matrix given by
2
3 20 0:25 19:75 6 7 H ¼ 4 20 20:25 0:25 5: 20 19:75 0:25 Choosing the function function, we get
y ¼ A0 þ
1 P
ð23Þ
sinðinxÞ expði1xÞ as the appropriate basis
i¼0
1 X Ai sinðinxÞ expði1xÞ;
ð24Þ
i¼1
where An is a vector in general, 1 and n are the constants, we take 1 = 1/4 and n = 1/2 in this example. When x = 0, the following set of equations can be derived
9 1 1 X X > > Ai ðinÞy00 ð0Þ ¼ 2Ai ði1ÞðinÞ > > > > = i¼1 i¼1 " # > ðnþ1Þ=2 1 > X X > n n!ð1Þjþ1 > Ai i 1n2jþ1 n2j1 > yðnÞ ð0Þ ¼ > ð2j1Þ!ðn2jþ1Þ! ;
repeating above steps, the approximate numerical solution of the given stiff system can be obtained. The numerical results are illustrated in Table 3.
4. Conclusions In this paper, according to the characteristic of stiff systems, an explicit method which applies an appropriate basis function to obtain the approximate numerical solution has been presented. The method is very simple and effective for most of differential equation systems, and it has been applied to the solution of stiff systems of point reactor neutron kinetics equations. Numerical examples show that the approach is promising and the research is worth to continue in this direction.
Acknowledgments This research is supported by the National Natural Science Foundation of China (Project No. 11301540) and the Natural Science Foundation of Naval University of Engineering.
yð0Þ ¼ A0 y0 ð0Þ ¼
i¼1
ð25Þ
j¼1
Also from Eq. (21) following equation can be obtained
yðnÞ ðxÞ ¼ Hyðn1Þ ðxÞ
ð26Þ
Applying the initial conditions yðxÞjx¼0 ¼ y0 , the values of y0 ð0Þ; y00 ð0Þ; ; yðnÞ ð0Þ can be obtained. Then according to Eq. (25), if the m-order approximation solution is chosen, the following system of equations will be obtained
" ðnþ1Þ=2 m X X n Ai i i¼1
j¼1
# n!ð1Þjþ1 n2jþ1 2j1 ¼ yðnÞ ð0Þ: 1 n ð2j 1Þðn 2j þ 1Þ!
ð27Þ
When n ¼ 1; 2; . . . ; m, according to Eq. (27), the linear equations of Ai can be gotten, and it is easy for us to develop a computer code to calculate the coefficient Ai. The corresponding mth-order approximation solution of the given stiff system is obtained as follows
y ¼ y0 þ
m X Ai sinðinxÞ expði1xÞ
ð28Þ
i¼1
Making step size of x to be h and substituting it into the mth-order approximation solution of y, y at x = x0 + h can be obtained. And
References Abdallah, A., Nahla, A., 2011. An efficient technique for the point reactor kinetics equations with Newtonian temperature feedback effects. Ann. Nucl. Energy 38, 2810–2817. Aboanber, A.E., 2004. On Padé approximations to the exponential function and application to the point kinetics equations. Prog. Nucl. Energy 44, 347–368. Aboanber, A.E., Hamada, Y.M., 2003. Power series solution (PWS) of nuclear reactor dynamics with Newtonian temperature feedback. Ann. Nucl. Energy 30, 1111– 1122. Butcher, J.C., 2000. Numerical methods for ordinary differential equations in the 20th century. J. Comput. Appl. Math. 125, 1–29. Chen, W.Z., Hao, J.L., Chen, L., Li, H.F., 2013. Solution of point reactor neutron kinetics equations with temperature feedback by singularly perturbed method. Sci. Technol. Nucl. Install., 1–6. Article ID 261327. Guzel, N., Bayram, M., 2005. On the numerical solution of stiff systems. Appl. Math. Comput. 170, 230–236. Hajmohammadi, M.R., Nourazar, S.S., 2014a. Conjugate forced convection heat transfer from a heated flat plate of finite thickness and temperature-dependent thermal conductivity. Heat Transfer Eng. 35, 863–874. Hajmohammadi, M.R., Nourazar, S.S., 2014b. On the solution of characteristic value problems arising in linear stability analysis; Semi analytical approach. Appl. Math. Comput. 239, 126–132. Hajmohammadi, M.R., Nourazar, S.S., 2014c. On the insertion of a thin gas layer in micro cylindrical Couette flows involving power-law liquids. Int. J. Heat Mass Transfer 75, 97–108. Hajmohammadi, M.R., Nourazar, S.S., Manesh, A.H., 2012. Semi-analytical treatments of conjugate heat transfer. J. Mech. Eng. Sci. 227, 492–503. Hamada, Y.M., 2013. Confirmation of accuracy of generalized power series method for the solution of point kinetics equations with feedback. Ann. Nucl. Energy 55, 184–193.
W. Chen et al. / Annals of Nuclear Energy 75 (2015) 353–357 Koclas, J., Sissaoui, M.T., Hébert, A., 1996. Solution of the improved and generalized quasi-static methods using an analytic calculation or a semi-implicit scheme to compute the precursor equations. Ann. Nucl. Energy 23, 1127–1142. Li, H.F., Chen, W.Z., Luo, L., et al., 2009. A new integral method for solving the point reactor neutron kinetics equations. Ann. Nucl. Energy 36, 427–432. Li, H.F., Shang, X.L., Chen, W.Z., 2010. An accurate solution of point kinetics equations of one-group delayed neutrons and an extraneous neutron source for step reactivity insertion. Chin. Sci. Bull. 55, 4116–4119. Liao, S., Tan, Y., 2007. A general approach to obtain series solutions of nonlinear differential equations. Stud. Appl. Math. 119, 297–354.
357
Tashakor, S., Jahanfarnia, G., Hashemi-Tilehnoee, M., 2010. Numerical solution of the point reactor kinetics equations with fuel burn-up and temperature feedback. Ann. Nucl. Energy 37, 265–269. Wazwaz, A.M., 2010. The combined Laplace transform–Adomian decomposition method for handling nonlinear Volterra integro–differential equations. Appl. Math. Comput. 216, 1304–1309. Zhu, Q., Shang, X.L., Chen, W.Z., 2012. Homotopy analysis solution of point reactor kinetics equations with six-group delayed neutrons. Acta Physica Sinica 61, 070201-1–070201-5.
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Explicit appropriate basis function method for numerical solution of stiff systems Wenzhen Chen a,b,⇑, Hongguang Xiao a, Haofeng Li a, Ling Chen a a b
Department of Nuclear Energy Science and Engineering, Naval University of Engineering, Wuhan 430033, China Institute of Thermal Science and Power Engineering, Naval University of Engineering, Wuhan 430033, China
a r t i c l e
i n f o
Article history: Received 11 April 2014 Received in revised form 7 July 2014 Accepted 20 August 2014
Keywords: Stiff systems Basis function Ordinary differential equations Explicit method
a b s t r a c t In this paper, an explicit numerical method, called the appropriate basis function method, is presented. The explicit appropriate basis function method differs from the power series method because it employs an appropriate basis function such as the exponential function, or periodic function, other than a polynomial, to obtain approximate numerical solutions. The method is successful and effective for the numerical solution of the first order ordinary differential equations. Two examples are presented to show the ability of the method for dealing with linear and nonlinear systems of differential equations. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction Numerous works have been focusing on the development of more advanced and efficient methods for solving ordinary differential equations (ODEs) (Butcher, 2000; Hajmohammadi et al., 2012, 2014a). For example, Wazwaz (2010) presented a new algorithm for solving in ODE’s of integro–differential type, and Hajmohammadi and Nourazar (2014b, c) presented a new algorithm for solving in ODE’s of eigen-value type based on semi-analytical method and pure-analytical method for solving in ODE’s of non-linear type, respectively. Implicit numerical methods are usually used for solving ordinary differential equations, in particular for some stiff problems, in which the implicitness will cause the method very involved compared to explicit methods (Li et al., 2009; Zhu et al., 2012). Indeed, explicit numerical methods require smaller step sizes in situations where implicit methods would not. It should be emphasized that explicit methods also show important advantages, some authors have developed many explicit methods for stiff systems (Guzel and Bayram, 2005). There has been a great deal of research that focused on eliminating the stiffness problem of reactor kinetics (Koclas et al., 1996; Chen et al., 2013). And there are several methods especially adapted for solving the initial value problems for stiff systems of ⇑ Corresponding author at: Department of Nuclear Energy Science and Engineering, Naval University of Engineering, Wuhan 430033, China. Tel.: +86 027 13871167436. E-mail address: [email protected] (W.Z. Chen). http://dx.doi.org/10.1016/j.anucene.2014.08.040 0306-4549/Ó 2014 Elsevier Ltd. All rights reserved.
ordinary differential equations (Aboanber and Hamada, 2003; Aboanber, 2004; Tashakor et al., 2010). Among the methods are numerical integration using Simpson’s rule, finite element method, Runge–Kutta procedures, quasi-static method, piecewise polynomial approach and other methods (Li et al., 2010; Abdallah and Nahla, 2011; Hamada, 2013). Most of these methods are successful in some specific problems, but still suffer, more or less, from disadvantages as mentioned by Chen et al., 2013. In this paper, an explicit numerical method for stiff systems is developed and tested. The method differs from the power series method (Guzel and Bayram, 2005) in the idea that employs an appropriate basis function, such as the exponential function, periodic function, and so on, other than a polynomial to obtain approximate numerical solutions. In the other hand, the method is considerably more accurate than the power series method, as demonstrated in following sections. 2. Appropriate basis function method 2.1. Basic ideas Since every ordinary differential equation with order n can be written as a system consisting of n ordinary differential equations with order one, our study is restricted to a system of first order differential equations, which can be written as follows
y0 ¼ f ðt; yÞ; t 2 ½0; T; kyk < þ1; yð0Þ ¼ y0 ;
ð1Þ
354
W. Chen et al. / Annals of Nuclear Energy 75 (2015) 353–357
whose theoretical solution is y(t). Let yn be an approximation to y(tn). Here
terms, Eq. (4) can also give the very accurate approximation. For example, when m = 2, Eq. (4) can be written as follows
y ¼ ½y1 ; y2 ; . . . ; yk T ;
tanhðtÞ 1 2e2t þ 2e4t e5t ;
T
f ¼ ½f 1 ; f 2 ; . . . ; f k ; and
y0 ¼ ½y01 ; y02 ; . . . ; y0k T : It is assumed that f and y are sufficient differentiability (Guzel and Bayram, 2005). The power series method (Guzel and Bayram, 2005) assumed that the solution of Eq. (1) can be expressed by a polynomial expression as follows
y¼
1 X Ai xi ;
ð2Þ
i¼0
where Ai is a vector function which is the same size as y0. However, the authors think that it should choose an appropriate function as the basis function according to the characteristic of stiff systems in order to improve the computational efficiency and accuracy. It is well known that a real function can be expressed by some different basis functions. For example, let us consider u(t) = tanh (t), which can be expressed by a polynomial expression as follows
tanhðtÞ t
m!þ1
1 X Ai xi ðxÞ;
ð5Þ
i¼0
where xi(x) is the ith term of the basis functions. For example, if the set of appropriate basis functions is fexpðixÞji 0g , xi(x) will be exp (ix).
ð3Þ
Then Eq. (3) converges to the exact solution u(t) only in a small region 0 6 t < 0.5 (see Fig. 1). However, by means of the exponential function as the basis function, one has
tanhðtÞ 1 þ lim 2
y¼
2.2. Scheme of the method
t 3 2t 5 17t 7 þ þ ; 3 15 315
"
which agrees well with the exact solution tanh (t) in the whole region 0 6 t < + 1. In addition, the numerical tests in Table 1 demonstrate that the efficiency of the exponential function is better than that of a polynomial expression. When the approximation solution is expressed by the exponential function u(t), it will take a computing time of 0.094 s to obtain the value of u(t) in the time interval [0, 3/2] with the time step h = 0.00001 s. However, using a polynomial as the basis function, it will take a computing time of 0.203 to 0.297 s. Therefore, one can get the best approximation by means of an appropriate basis function. For example, a periodic solution is expressed more efficiently by periodic basis functions than by a polynomial expression. The solution of Eq. (1) can be expressed by a set of appropriate basis functions fxi ðiÞji 0g as follows
#
m X
ð1Þn e2nt þ ð1Þmþ1 eð2mþ1Þt ;
The way that the appropriate basis function method is used in practice is carried out by the following computation steps: Step 1: Choosing the order of approximation solution, m. Step 2: When x = 0, according to Eq. (5), we can get
ð4Þ
yðnÞ ð0Þ ¼
n¼1
which converges to the exact solution u(t) in the whole region 0 6 t < + 1 (Liao and Tan, 2007). And, even taking the first few
m X ðnÞ Ai xi ð0Þ:
ð6Þ
i¼0
Step 3: According to Eq. (1) and the initial conditions, the values or expressions of yð0Þ; y0 ð0Þ; y00 ð0Þ; . . . ; yðmÞ ð0Þ can be obtained by the iterative computation. Step 4: When n ¼ 0; 1; 2; . . . ; m, according to Eq. (6), the linear equations of Ai can be derived, and it is easy for us to develop a computer code to calculate the coefficient Ai. Then, substituting Ai into Eq. (5), the m th-order approximation solution of Eq. (1) is obtained . Step 5: Making step size of x to be h and substituting it into the mth-order approximation solution of Eq. (1), we have y at x = x0 + h. Step 6: Repeating above steps 1 to 5, we can obtain the numerical solution of Eq. (1). 3. Numerical experiments The appropriate basis function method and the power series method (Guzel and Bayram, 2005) are compared in the following cases. We arranged the independent variable x, the step size h,
Fig. 1. Exact and approximation solutions of u(t) = tanh (t).
Table 1 Comparison of approximation solutions of u(t) = tanh (t) in the time interval [0, 3/2] by using h = 105. Type of approximation
Solutions for different time t = 0.8
t=1
Computing time (s) t = 1.1
t = 1.2
t = 1.3
t = 1.4
t = 1.5
uðtÞ ¼ t t3 þ 2t5 17t 315
3
5
0.7491
1.0127
1.1954
1.4260
1.7142
2.0677
2.4904
0.297
t3 3
5 þ 2t5 2t
0.7604
1.0667
1.3005
1.6193
2.0528
2.6366
3.4125
0.203
0.6594 t = 0.8 0.6640
0.7592 t=1 0.7616
0.7989 t = 1.1 0.8005
0.8326 t = 1.2 0.8337
0.8610 t = 1.3 0.8617
0.8849 t = 1.4 0.8854
0.9048 t = 1.5 0.9052
0.094
uðtÞ ¼ t
7
þ 2e4t e5t uðtÞ ¼ 1 2e Accurate solution u(t) = tanh (t)
355
W. Chen et al. / Annals of Nuclear Energy 75 (2015) 353–357 Table 2 Comparison of theoretical and numerical results in problem 1. Numerical results h
N
y
Theoretical results
Appropriate basis function method (i = 4)
1 ¼ 1 0.05
1
1
1
y1 y2 y1 y2
0.90483741803596 0.95122942450071 0.13533528323661 0.36787944117144
Power series method (Guzel and Bayram, 2005) (i = 4)
1 ¼ 0:6
Solution
Error
Solution
Error
Solution
Error
0.90483741803596 0.95122942450071 0.13533528323661 0.36787944117144
0 0 0 0
0.90483743957714 0.95122944625838 0.13611513643068 0.36883051686970
2.1541e8 2.1758e8 7.7985e4 9.5108e4
0.90483750000000 0.95122942708333 0.33333333333333 0.37500000000000
8.1964e8 2.5826e9 0.197998 0.007121
the computation step number N, the theoretical solution, and the numerical results by the appropriate basis function method and the power series method, respectively. All computations are carried out by the MATLAB code with double precision arithmetic.
y001 ð0Þ ¼ 1002y01 ð0Þ þ 2000y2 ð0Þy02 ð0Þ ¼ 4 y002 ð0Þ ¼ y01 ð0Þ ð1 þ 2y2 ð0ÞÞy02 ð0Þ ¼ 1
00 02 00 y000 1 ð0Þ ¼ 1002y1 ð0Þ þ 2000y2 ð0Þ þ 2000y2 ð0Þy2 ð0Þ ¼ 8 00 00 02 y000 2 ð0Þ ¼ y1 ð0Þ ð1 þ 2y2 ð0ÞÞy2 ð0Þ 2y2 ð0Þ ¼ 1
3.1. Problem 1
ð16Þ
;
:
ð17Þ
Then, the linear equation can be given as follows The test problem considers the following system of differential equations
where
y01 ¼ 1002y1 þ 1000y22
with initial conditions: y1(0) = 1 and y2(0) = 1. The exact solutions are y1(x) = e2x and y2(x) = ex. y1 and y2 can be represented by the exponential basis function other than a polynomial as follows, respectively 1 X
a1;i e1ix and y2 ðxÞ ¼
i¼1
1 X a2;i e1ix ;
ð8Þ
i¼1
where 1 is a constant, it can provide us with great freedom to choose appropriate basis functions to approximate the stiff system. Then, the n th derivatives of y1 and y2 are given by 1 1 X X n n ðnÞ ðnÞ y1 ðxÞ ¼ a1;i ð1iÞ e1ix and y2 ðxÞ ¼ a2;i ð1iÞ e1ix : i¼1
ð9Þ
i¼1
When m = 4, the 4th-order approximation solutions of Eq. (7) are given by
y1 ðxÞ ¼
4 X
a1;i e1ix
4 X and y2 ðxÞ ¼ a2;i e1ix ;
i¼1
2
ð7Þ
y02 ¼ y1 y2 ð1 þ y2 Þ;
y1 ðxÞ ¼
Aa ¼ B;
6 60 6 6 61 6 6 60 6 A¼6 6 12 6 6 60 6 6 6 13 4 0
i¼1
ðnÞ
i¼1
y2 ð0Þ ¼ 1
y02 ð0Þ
¼ y1 ð0Þ y2 ð0Þð1 þ y2 ð0ÞÞ ¼ 1
0
0
1
1
1
1
1
1
0
0
0
0
0
0
1
1
1
2
2
2
0
0
0
0
2
1
2
1
12
13 13 13
0
0
0
0
13 13 13 13
1
0
0
0
1
0
3
7 17 7 7 07 7 7 17 7 7; 07 7 7 27 1 7 7 07 5
19 2
1 2
ð18Þ
1 41 1 ; 4 þ 71 1
2
3
1
Finally, the 4th-order approximation solutions are obtained. Note that, when 1 ¼ 1, the exact solutions y1(x) = e2x and y2(x) = ex can also be obtained. The numerical results are illustrated in Table 2. 3.2. Problem 2
ð13Þ
Let us consider the following system of differential equations
i¼1
y01 ð0Þ ¼ 1002y1 ð0Þ þ 1000y22 ð0Þ ¼ 2
0
1
0
þ 72 12 þ 12 13 ; 1 11 11 12 16 13 T : 6
ð11Þ
ð12Þ
ð14Þ
;
0
1612 413 ;4 þ 1411 þ 1412 þ 413 ;
6
And according to Eq. (7) and the initial conditions, following equations can be derived
y1 ð0Þ ¼ 1
0
a ¼ 4 þ 26 11 þ 612 þ 43 13 ; 6 1911 3
i¼1
4 4 X X n n ðnÞ a1;i ð1iÞ and y2 ð0Þ ¼ a2;i ð1iÞ ;
0
11 412 43 13 ; 4 þ 133 11 þ 23 12 þ 16 13 ; 1 11 3
y1 ð0Þ ¼ a1;1 þ a1;2 þ a1;3 þ a1;4 and y2 ð0Þ
y1 ð0Þ ¼
1
Solving this linear equation results in
ð10Þ
Substituting x = 0 into Eqs. (10) and (11) results in
¼ a2;1 þ a2;2 þ a2;3 þ a2;4 ;
1
a ¼ ½a1;1 ; a1;2 ; a1;3 ; a1;4 ; a2;1 ; a2;2 ; a2;3 ; a2;4 T :
i¼1
4 4 X X n n ðnÞ ¼ a1;i ð1iÞ e1ix and y2 ðxÞ ¼ a2;i ð1iÞ e1ix ;
1
B ¼ ½1; 1; 2; 1; 4; 1; 8; 1T ;
According Eq. (10), following equations can be gotten ðnÞ y1 ðxÞ
1
y01 y02 y03
¼ 20y1 0:25y2 19:75y3 ; y1 ð0Þ ¼ 1; ¼ 20y1 20:25y2 þ 0:25y3 ; ¼ 20y1 19:75y2 0:25y3 ;
y2 ð0Þ ¼ 0; y3 ð0Þ ¼ 1:
ð19Þ
The analytical solutions of Eq. (19) are
y1 ¼ ½e1=2x þ e20x ðcosð20xÞ þ sinð20xÞÞ=2; ;
ð15Þ
y2 ¼ ½e1=2x e20x ðcosð20xÞ sinð20xÞÞ=2; y3 ¼ ½e1=2x þ e20x ðcosð20xÞ sinð20xÞÞ=2:
ð20Þ
356
W. Chen et al. / Annals of Nuclear Energy 75 (2015) 353–357
Table 3 Comparison of theoretical and numerical results in problem 2. Numerical results h
N
0.01
20
0.1
2
0.025
y
y1 y2 y3 y1 y2 y3 y1 y2 y3 y1 y2 y3
20
0.1
Theoretical results
5
0.43950209815007 0.45147399867173 0.45336341936422 0.43950209815007 0.45147399867173 0.45336341936422 0.38936899538135 0.38940708916983 0.38939369390157 0.38936899538135 0.38940708916983 0.38939369390157
Appropriate basis function method (m = 6)
Power series method (Guzel and Bayram, 2005) (i = 6)
Solution
Error
Solution
Error
0.43950209807394 0.45147411086509 0.45336330759892 0.45326849690173 0.44691922771744 0.45791817689314 0.38936955757924 0.38940855772086 0.38939222627783 0.38939761019105 0.38939766195830 0.38940310409674
7.61273e011 1.12193e007 1.11765e007 0.013766398 0.004554770 0.004554757 5.62197e007 1.46855e006 1.46762e006 2.86148e005 9.42721e006 9.41019e006
0.43950209127250 0.45147399307674 0.45336342495922 0.42698661025269 0.51933228926504 0.38550512877121 0.38936730990046 0.38940763611512 0.38939314695629 0.39162253788905 0.38878523056801 0.39001555250403
6.87757e009 5.59499e009 5.59499e009 0.012515487 0.067858290 0.067858290 1.68548e006 5.46945e007 5.46945e007 0.002253542 6.21858e004 6.21858e004
In addition, Eq. (19) can be rewritten by a matrix formulation as follows
dy ¼ Hy; dx
ð21Þ
yðxÞjx¼0 ¼ y0 ;
ð22Þ
where y = [y1, y2, y3]T and y0 = [1, 0, 1]T. H is a square matrix given by
2
3 20 0:25 19:75 6 7 H ¼ 4 20 20:25 0:25 5: 20 19:75 0:25 Choosing the function function, we get
y ¼ A0 þ
1 P
ð23Þ
sinðinxÞ expði1xÞ as the appropriate basis
i¼0
1 X Ai sinðinxÞ expði1xÞ;
ð24Þ
i¼1
where An is a vector in general, 1 and n are the constants, we take 1 = 1/4 and n = 1/2 in this example. When x = 0, the following set of equations can be derived
9 1 1 X X > > Ai ðinÞy00 ð0Þ ¼ 2Ai ði1ÞðinÞ > > > > = i¼1 i¼1 " # > ðnþ1Þ=2 1 > X X > n n!ð1Þjþ1 > Ai i 1n2jþ1 n2j1 > yðnÞ ð0Þ ¼ > ð2j1Þ!ðn2jþ1Þ! ;
repeating above steps, the approximate numerical solution of the given stiff system can be obtained. The numerical results are illustrated in Table 3.
4. Conclusions In this paper, according to the characteristic of stiff systems, an explicit method which applies an appropriate basis function to obtain the approximate numerical solution has been presented. The method is very simple and effective for most of differential equation systems, and it has been applied to the solution of stiff systems of point reactor neutron kinetics equations. Numerical examples show that the approach is promising and the research is worth to continue in this direction.
Acknowledgments This research is supported by the National Natural Science Foundation of China (Project No. 11301540) and the Natural Science Foundation of Naval University of Engineering.
yð0Þ ¼ A0 y0 ð0Þ ¼
i¼1
ð25Þ
j¼1
Also from Eq. (21) following equation can be obtained
yðnÞ ðxÞ ¼ Hyðn1Þ ðxÞ
ð26Þ
Applying the initial conditions yðxÞjx¼0 ¼ y0 , the values of y0 ð0Þ; y00 ð0Þ; ; yðnÞ ð0Þ can be obtained. Then according to Eq. (25), if the m-order approximation solution is chosen, the following system of equations will be obtained
" ðnþ1Þ=2 m X X n Ai i i¼1
j¼1
# n!ð1Þjþ1 n2jþ1 2j1 ¼ yðnÞ ð0Þ: 1 n ð2j 1Þðn 2j þ 1Þ!
ð27Þ
When n ¼ 1; 2; . . . ; m, according to Eq. (27), the linear equations of Ai can be gotten, and it is easy for us to develop a computer code to calculate the coefficient Ai. The corresponding mth-order approximation solution of the given stiff system is obtained as follows
y ¼ y0 þ
m X Ai sinðinxÞ expði1xÞ
ð28Þ
i¼1
Making step size of x to be h and substituting it into the mth-order approximation solution of y, y at x = x0 + h can be obtained. And
References Abdallah, A., Nahla, A., 2011. An efficient technique for the point reactor kinetics equations with Newtonian temperature feedback effects. Ann. Nucl. Energy 38, 2810–2817. Aboanber, A.E., 2004. On Padé approximations to the exponential function and application to the point kinetics equations. Prog. Nucl. Energy 44, 347–368. Aboanber, A.E., Hamada, Y.M., 2003. Power series solution (PWS) of nuclear reactor dynamics with Newtonian temperature feedback. Ann. Nucl. Energy 30, 1111– 1122. Butcher, J.C., 2000. Numerical methods for ordinary differential equations in the 20th century. J. Comput. Appl. Math. 125, 1–29. Chen, W.Z., Hao, J.L., Chen, L., Li, H.F., 2013. Solution of point reactor neutron kinetics equations with temperature feedback by singularly perturbed method. Sci. Technol. Nucl. Install., 1–6. Article ID 261327. Guzel, N., Bayram, M., 2005. On the numerical solution of stiff systems. Appl. Math. Comput. 170, 230–236. Hajmohammadi, M.R., Nourazar, S.S., 2014a. Conjugate forced convection heat transfer from a heated flat plate of finite thickness and temperature-dependent thermal conductivity. Heat Transfer Eng. 35, 863–874. Hajmohammadi, M.R., Nourazar, S.S., 2014b. On the solution of characteristic value problems arising in linear stability analysis; Semi analytical approach. Appl. Math. Comput. 239, 126–132. Hajmohammadi, M.R., Nourazar, S.S., 2014c. On the insertion of a thin gas layer in micro cylindrical Couette flows involving power-law liquids. Int. J. Heat Mass Transfer 75, 97–108. Hajmohammadi, M.R., Nourazar, S.S., Manesh, A.H., 2012. Semi-analytical treatments of conjugate heat transfer. J. Mech. Eng. Sci. 227, 492–503. Hamada, Y.M., 2013. Confirmation of accuracy of generalized power series method for the solution of point kinetics equations with feedback. Ann. Nucl. Energy 55, 184–193.
W. Chen et al. / Annals of Nuclear Energy 75 (2015) 353–357 Koclas, J., Sissaoui, M.T., Hébert, A., 1996. Solution of the improved and generalized quasi-static methods using an analytic calculation or a semi-implicit scheme to compute the precursor equations. Ann. Nucl. Energy 23, 1127–1142. Li, H.F., Chen, W.Z., Luo, L., et al., 2009. A new integral method for solving the point reactor neutron kinetics equations. Ann. Nucl. Energy 36, 427–432. Li, H.F., Shang, X.L., Chen, W.Z., 2010. An accurate solution of point kinetics equations of one-group delayed neutrons and an extraneous neutron source for step reactivity insertion. Chin. Sci. Bull. 55, 4116–4119. Liao, S., Tan, Y., 2007. A general approach to obtain series solutions of nonlinear differential equations. Stud. Appl. Math. 119, 297–354.
357
Tashakor, S., Jahanfarnia, G., Hashemi-Tilehnoee, M., 2010. Numerical solution of the point reactor kinetics equations with fuel burn-up and temperature feedback. Ann. Nucl. Energy 37, 265–269. Wazwaz, A.M., 2010. The combined Laplace transform–Adomian decomposition method for handling nonlinear Volterra integro–differential equations. Appl. Math. Comput. 216, 1304–1309. Zhu, Q., Shang, X.L., Chen, W.Z., 2012. Homotopy analysis solution of point reactor kinetics equations with six-group delayed neutrons. Acta Physica Sinica 61, 070201-1–070201-5.