- Email: [email protected]
A non-standard explicit integration scheme for initial-value problems
Applied Mathematics and Computation 189 (2007) 710–718 www.elsevier.com/locate/amc
A non-standard explicit integration scheme for initial-value problems Higinio Ramos Scientific Computing Group, Universidad de Salamanca, Spain Escuela Polite´cnica Superior de Zamora, Campus Viriato, 49022 Zamora, Spain
Abstract In this paper we present the construction of a non-standard explicit algorithm for initial-value problems. The method results to be of second order and A-stable. This new algorithm has been proven to be suitable for solving different kind of initial-value problems, specifically, non-singular problems, singular problems, stiff problems and singularly perturbed problems. Some numerical experiments are considered in order to check the behaviour of the method when applied to a variety of initial-value problems. Ó 2006 Elsevier Inc. All rights reserved. Keywords: Non-standard algorithm; Stiff initial-value problems; Singular initial-value problems; Singularly perturbed initial-value problems
1. Introduction A non-standard numerical integration method for initial-value problems (IVP) is presented. It is shown that the new method has second-order accuracy and is A-stable. This favorable property of A-stability suggests that the new scheme may be used for solving stiff systems. Moreover, the non-linear characteristic of the method suggests that it may also be suitable to solve special problems for which conventional numerical integrators result inefficient, such as singular IVPs or singularly-perturbed IVPs. A singular IVP is characterized by the presence of a pole in the solution or by a discontinuous low order derivative of the solution. A simple example is the innocent looking IVP y0 ¼ 1 þ y2;
yð0Þ ¼ 1;
x 2 ½0; b;
ð1Þ
with b > p=4. The theoretical solution is yðxÞ ¼ tanðx þ p=4Þ which has a simple pole at x ¼ p=4. Singularly-perturbed IVPs usually contain a small parameter that multiplies the first-order derivative. These problems are characterized by the presence of thin layers where the solution varies very fast, whereas away from the layers the solution behaves regularly and varies slowly. For the analytical and numerical treatment of this E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.11.134
H. Ramos / Applied Mathematics and Computation 189 (2007) 710–718
711
type of problems the reader is referred to the books by O’Malley [1] and Roos et al. [2]. A common strategy for dealing with this type of problems consists of dividing the domain of integration in two subdomains and then to apply a different scheme on each subdomain (see [3,4]). The use of classical explicit numerical methods for these type of problems requires very small time steps. And the use of standard implicit methods may require a large number of iterations for convergence within each time step. The method presented in this article has been tested on a variety of IVPs and the results show that it may be appropriate for solving IVPs of different nature. The paper has been organized as follows. In Section 2 the non-standard algorithm is developed, in Sections 3 and 4 the main characteristics of the method are established. A section of numerical results on different problems is presented to show the performance of the new method. And, finally, a summary of the main conclusions puts an end to the paper. 2. Non-standard explicit method Consider the IVP y 0 ¼ f ðx; yÞ; yðaÞ ¼ y 0 ; y; f ðx; yÞ 2 R; x 2 ½a; b R;
ð2Þ
where f is assumed to satisfy all the requirements in order that (2) has a unique solution. The interval ½a; b is divided into a number of subintervals ½xj ; xjþ1 with x0 ¼ a and xj ¼ x0 þ jh, being h the step size. Suppose we have solved numerically the problem in (2) up to a point xn, and have obtained a value yn as an approximation of y(xn). Assuming the localization hypothesis [6], yðxn Þ ¼ y n , we are interested in obtaining an approximate value, yn+1, for the true value yðxnþ1 Þ. For that purpose, following the ideas in [5], we suggest an approximation to the theoretical solution yðxn þ hÞ of (2) given by yðxn þ hÞ ¼ yðxn Þ þ
hy 0 ðxn Þ ; aðhÞ þ yðxn Þ
ð3Þ
where a(h) is a sufficiently differentiable unknown function of the step size that has to be determined and it is assumed that aðhÞ þ yðxn Þ 6¼ 0. From (3) it results that F n ðy; a; hÞ ¼ ðyðxn þ hÞ yðxn ÞÞðaðhÞ þ yðxn ÞÞ hy 0 ðxn Þ ¼ 0: Now, if we consider a(h) expanded in Taylor series about h ¼ 0, after expanding F n ðy; a; hÞ in Taylor series about xn the following expression is obtained: 1 F n ðy; a; hÞ ¼ ½1 þ að0Þ þ yðxn Þy 0 ðxn Þh þ ½2a0 ð0Þy 0 ðxn Þ þ ðað0Þ þ yðxn ÞÞy 00 ðxn Þh2 þ Oðh3 Þ: 2
ð4Þ
Imposing that the coefficients of h and h2 in (4) vanish, we obtain a system of equations from which it is readily deduced that að0Þ ¼ 1 yðxn Þ; a0 ð0Þ ¼
y 00 ðxn Þ ; 2y 0 ðxn Þ
provided that y 0 ðxn Þ 6¼ 0. Introducing the above values in the Taylor series of a(h) we have aðhÞ ¼ 1 yðxn Þ
y 00 ðxn Þ h þ Oðh2 Þ: 2y 0 ðxn Þ
ð5Þ
From (3) and (5) it is readily deduced the numerical scheme, which may be written in the form y nþ1 ¼ y n þ
2hf 2n ; 2f n hf 0n
ð6Þ
712
H. Ramos / Applied Mathematics and Computation 189 (2007) 710–718
where y n ¼ yðxn Þ; y nþ1 ’ yðxnþ1 Þ; fn ¼ f ðxn ; y n Þ and fn0 ¼
of of ðxn ; y n Þ þ ðxn ; y n Þfn : ox oy
In particular, when the differential equation in (2) is autonomous, y 0 ¼ f ðyÞ, the method in (6) reads y nþ1 ¼ y n þ
2hfn ; 2 hjn
ð7Þ
ðy n Þ. where we use now the notation fn ¼ f ðy n Þ and jn ¼ of oy Thus we have developed an explicit non-linear one-step method of the form y nþ1 ¼ y n þ hUf ðy n ; xn ; hÞ, where Uf ðy n ; xn ; hÞ is the incremental function, and the subscript f on the right hand side indicates that the dependence of U on yn is through the function f (see [6]). After applying this scheme we will take as approximation for the true solution of (2) at xn+1 the value yn+1 given respectively by (7) or either (6) depending on whether the problem is autonomous or not. Repeating the procedure along the nodes on the integration interval we will obtain a discrete solution for the problem in (2). Non-linear multistep methods are usually designed for dealing with unconventional problems such as stiff or singular initial-value problems, for which the classical schemes generally perform poorly [7,9–12]. On the contrary, some of the non-linear schemes specifically designed for singular problems do not perform well on non-singular problems [12]. In the numerical section we will see that the proposed non-linear scheme is suitable for a variety of IVPs of different nature. 3. Local truncation error It follows from the construction of the method in (6) that the new scheme is at least of second order. In fact, if we consider the functional given by LðzðxÞ; hÞ ¼ zðx þ hÞ zðxÞ
2hz0 ðxÞ2 ; 2z0 ðxÞ hz00 ðxÞ
ð8Þ
where z(x) is an arbitrary function defined on ½a; b and differentiable as often as we need, after expanding in Taylor series about x and collecting terms in h we obtain ! 2 z00 ðxÞ zð3Þ ðxÞ 3 LðzðxÞ; hÞ ¼ þ ð9Þ h þ Oðh4 Þ; 6 4z0 ðxÞ which means that the method has at least second order of accuracy. The local truncation error of the method may be written as ! 2 ðy 00n Þ y 000 n T nþ1 ¼ þ h3 þ Oðh4 Þ; 4y 0n 6 where y 0n ; y 00n ; y 000 n denote, respectively the numerical approximations to the first, second and third derivatives of y(x) at the point xn. Note that if we solve the resulting equation from equating to zero the coefficient of h3 in (9) it may be verified easily that the method in (6) is exact when the solution of the differential equation in (2) is of the form yðxÞ ¼ a þ
b ; cþx
where a; b; c are real constants. This occurs for example for each of the following differential equations: 2
ðaÞ
2y 0 ðxÞy ð3Þ ðxÞ ¼ 3y 00 ðxÞ ;
ðbÞ
2y 0 ðxÞ þ ðc þ xÞy 00 ðxÞ ¼ 0;
ðcÞ
4y 0 ðxÞ þ by 00 ðxÞ ¼ 0;
3
2
H. Ramos / Applied Mathematics and Computation 189 (2007) 710–718
713
2
ðdÞ 2y 0 ðxÞ þ ay 00 ðxÞ yðxÞy 00 ðxÞ ¼ 0; ðeÞ b þ ðc2 þ 2cx þ x2 Þy 0 ðxÞ ¼ 0; ðfÞ yðxÞ þ ðc þ xÞy 0 ðxÞ ¼ a; ðgÞ a2 2ayðxÞ þ yðxÞ2 þ by 0 ðxÞ ¼ 0:
ð10Þ
Remark 1. We observe that when the solution of the differential equation in (2) verifies the equation 2
2y 0 ðxÞ ¼ yðxÞy 00 ðxÞ; then we obtain as a particular case of our method the scheme given by y nþ1 ¼
y 2n ; y n hf ðy n Þ
that is, the first-order scheme in [8], which results to be exact when the solution of the differential equation is of a the form yðxÞ ¼ xþb ; a; b 2 R. This occurs for example in the scheme in (7) for the autonomous problem in the case that y n jn ¼ 2f ðy n Þ holds. 4. Stability analysis The stability of the scheme in (6) is examined by applying it to the Dahlquist’s test problem y 0 ¼ ky; ReðkÞ < 0, yielding the difference equation y nþ1 ¼
2 þ kh y : 2 kh n
Setting z ¼ kh the stability function is obtained as RðzÞ ¼
2þz ; 2z
which is the (1, 1)-Pade´ approximation to the exponential ez, and thus the method is A-stable (see [13]). The stability region of the method consists of the left-half complex plane. 5. Numerical results To test the proposed method, we have integrated several initial value-problems of different nature. All the numerical experiments were performed in double precision using a FORTRAN compiler on a personal computer with a Pentium based processor (2019 MHz). 5.1. A Riccati equation The first example we have considered consists of the Riccati equation coupled with the initial value given by 2
y 0 ðxÞ ¼ 4 þ 4yðxÞ yðxÞ ;
yð0Þ ¼ 1;
whose exact solution is yðxÞ ¼
2x 1 ; x1
and has a pole in x ¼ 1. Although this is a severe difficult problem for most of the usual integrators our method is capable of solving the problem exactly, that is to say, without local truncation error. We have integrated this problem with our method on the interval [0, 2] taking different values for the constant stepsize.
714
H. Ramos / Applied Mathematics and Computation 189 (2007) 710–718
In Table 1 the results for this problem are shown. The errors have been defined as the maximum of the absolute errors on the nodal points in the integration interval Emax ¼ max fjyðxj Þ y j jg; xj 2½a;b
or as the absolute error at the final point Eðx ¼ bÞ ¼ jyðbÞ y NI j; where NI refers to the number of steps used in the integration. The CPU column refers to the computation time measured in milliseconds. We observe that the method has no difficulties in crossing the point x ¼ 1. A plot of the absolute errors is shown in Fig. 1. The differential equation in this problem coincides with the type (g) in (10), and so the errors are only due to the effect of rounding errors. We observe that the bigger the number of steps is so it is the corresponding error. Remark 2. It is possible that for certain combinations of the values appearing in the denominators in (6) or either in (7), some of these denominators vanish. In that case, it is necessary to modify slightly the step size on this step and take mh instead of h, where the factor m is taken nearly unity, say m ¼ 0:95. 5.2. A stiff equation Stiff differential equations involve rapidly increasing or decaying transient solution, which results in a difficulty for most of the numerical integrators. The method in (6) was tested on the stiff problem taken from [14] y 0 ðxÞ ¼ 100yðxÞ þ 99e2x ;
yð0Þ ¼ 0;
which has exact solution yðxÞ ¼
33 2x ðe e100x Þ: 34
The problem has been integrated on the interval [0, 0.5] and the results are presented in Table 2 for different number of steps, NI. Table 1 Errors for y 0 ðxÞ ¼ 4 þ 4yðxÞ yðxÞ2 ; yð0Þ ¼ 1 NI 16 32 64 128 256 512
Eðx ¼ 2Þ
Emax 1:42 10 5:18 1013 2:70 1013 1:99 1011 1:55 1010 1:65 1011
0.00 8:88 1016 3:10 1015 3:55 1015 5:32 1015 2:79 1014
2.5·10
-13
2·10
-13
1.5·10
-13
1·10
-13
5·10
-14
CPU 14
0.5
1
1.5
Fig. 1. Absolute errors for Problem 5.1.
000 000 000 000 000 000
2
H. Ramos / Applied Mathematics and Computation 189 (2007) 710–718
715
Table 2 Errors for y 0 ðxÞ ¼ 100yðxÞ þ 99e2x ; yð0Þ ¼ 0 NI
Eðx ¼ 0:5Þ
Emax
CPU
16 32 64 128 256 512
3:12 102 2:06 105 1:74 106 3:28 107 7:30 108 1:73 108
2:31 101 7:81 102 1:78 102 4:14 103 1:03 103 2:55 104
000 000 000 000 000 016
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5 0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
Fig. 2. Exact and numerical solutions (NI ¼ 64) for Problem 5.2.
In Fig. 2 the exact and the numerical solutions are shown. For NI ¼ 64 we do not appreciate notably differences on the two plots. 5.3. A stiff system The above method may be also applied to a system of equations. If we have y; f ðx; yÞ 2 Rm in (2) we have just to consider the formula in (6) for every component. Let be the stiff system taken from [15] y 01 ðxÞ ¼ 1002y 1 ðxÞ þ 1000y 2 ðxÞ2 ; y 02 ðxÞ
¼ y 1 ðxÞ y 2 ðxÞð1 þ y 2 ðxÞÞ;
y 1 ð0Þ ¼ 1; y 2 ð0Þ ¼ 1;
which will be integrated with our method on the interval [0, 1]. The exact solution is y 1 ðxÞ ¼ e2x ;
y 2 ðxÞ ¼ ex :
The results for every component of the solution appear in Table 3.
Table 3 Errors for Problem 5.3 NI
Eðx ¼ 1Þ y 1 ðxÞ
y 2 ðxÞ
y 1 ðxÞ
y 2 ðxÞ
128 256 512 1024 2048
2:03 103 3:33 105 9:77 108 2:17 108 5:41 109
3:83 103 3:96 104 1:17 107 2:93 108 7:33 109
1:29 102 1:10 102 1:54 107 2:97 108 7:40 109
4:08 103 4:18 104 1:17 107 2:93 108 7:33 109
Emax
CPU
000 000 000 000 000
716
H. Ramos / Applied Mathematics and Computation 189 (2007) 710–718
5.4. A singular problem The next example we consider is the IVP given in (1). The solution becomes unbounded in the neighborhood of the singularity at x ¼ p=4 ’ 0:785398163397448. Finite difference methods based on local polynomial interpolation behave poorly if the solution has singularities. This problem has been used as a numerical test for different non-linear integrators intended for numerically solving singular initial-value problems [7–12]. In Table 4 we present the results for our method where we may observe the ability of the method to cross the singularity. In Fig. 3 the numerical solution after joining the points ðxn ; y n Þ is shown for NI ¼ 128. We observe that the numerical solution reproduces the behaviour of the true solution faithfully. 5.5. A singularly-perturbed problem Now we consider a singularly perturbed IVP given by 2
y 0 ðxÞ ¼ 2kyðxÞ ;
yð0Þ ¼ 10;
which has the exact solution yðxÞ ¼
10 : 1 þ 20kx
We have considered k ¼ 500 as in [16] so that the solution drops quickly from its initial value of 10 to very small values. This problem exhibits an initial layer near x ¼ 0. Table 5 shows the errors obtained with our non-linear method when the integration is performed on the interval [0, 1]. We note that even near x ¼ 0, the new method performs very well, even with the presence of the initial layer there. In fact, the differential equation in this problem coincides with the type (g) in (10), and thus the errors of the method are due only to round off considerations. In Fig. 4 the absolute errors are shown in logarithmic scale when NI ¼ 16.
Table 4 Errors for y 0 ðxÞ ¼ 1 þ yðxÞ2 ; yð0Þ ¼ 1 NI
Eðx ¼ 1Þ
Emax
CPU
32 64 128 256 512 1024 2048
7:18 103 1:79 103 4:48 104 1:12 104 2:80 105 7:00 106 1:75 106
13.9181 3.6385 1.2008 67.1306 16.9897 4.26056 1.06597
000 000 000 000 000 000 016
20 10
0.2
0.4
0.6
0.8
-10 -20 Fig. 3. Numerical solution (NI ¼ 128) for Problem 5.4.
1
H. Ramos / Applied Mathematics and Computation 189 (2007) 710–718
717
Table 5 Errors for y 0 ðxÞ ¼ 2kyðxÞ2 ; yð0Þ ¼ 10; k ¼ 500 NI
Eðx ¼ 1Þ
Emax
CPU
16 32 64 128
2:16 1019 6:50 1019 2:16 1019 0.00
9:71 1017 8:67 1016 7:63 1016 7:21 1016
000 000 000 016
Log (Err) -16 -16.5 -17 -17.5 -18 -18.5 0.2
0.4
0.6
0.8
1
x
Fig. 4. Absolute errors (NI ¼ 16) for Problem 5.5.
5.6. An autonomous problem The last problem to be considered is 3
y 0 ðxÞ ¼ 1 þ yðxÞ ;
yð0Þ ¼ 0;
which will be integrated using different stepsizes in the interval [0, 1.2]. The analytical solution is given implicitly by pffiffiffi pffiffiffi 2y 1 pffiffiffi þ 6 logð1 þ yÞ ¼ 3ð6x þ logð1 y þ y 2 ÞÞ; 3p þ 6 3 arctan 3 from which it is obtained using the Newton method that for x ¼ 1:2 it is yð1:2Þ ¼ 7:368587110472472161317. This value will be used as a reference for calculating the errors. In Table 6 the results obtained with the new method for this problem are presented. In Fig. 5 the numerical solution after joining the points ðxn ; y n Þ is shown for NI ¼ 256.
Table 6 Errors for y 0 ðxÞ ¼ 1 þ yðxÞ3 ; yð0Þ ¼ 0 NI
Eðx ¼ 1:2Þ
CPU
256 512 1024 2048 4096
5:8098 102 1:5670 102 4:0356 103 1:0107 103 2:4024 104
000 000 000 000 016
718
H. Ramos / Applied Mathematics and Computation 189 (2007) 710–718
7 6 5 4 3 2 1 0.2
0.4
0.6
0.8
1
1.2
Fig. 5. Numerical solution (NI ¼ 256) for Problem 5.6.
6. Conclusions In this article the development of a new non-linear explicit scheme for initial value problems is considered. The method is shown to have second order and to be A-stable. These characteristics suggested their use for solving singular and stiff initial-value problems. Numerical results demonstrate that the method can handle these kind of problems at low cost as it is shown by the CPU values in the tables. Acknowledgements The author wish to thank JCYL project SA 024/04 and CYT project BMF 2004-295 for financial support. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]
R.E. O’Malley, Singular Perturbation Methods for Ordinary Differential Equations, Springer-Verlag, New York, 1991. H.G. Roos, M. Stynes, L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations, Springer, Berlin, 1996. S.M. Roberts, A boundary-value technique for singular perturbation problems, J. Math. Anal. Appl. 79 (1982) 489–508. S. Natesan, J. Vigo-Aguiar, N. Ramanujam, A numerical algorithm for singular perturbation problems exhibiting weak boundary layers, Comput. Math. Appl. 45 (2003) 469–479. F.D. van Niekerk, Nonlinear one-step methods for initial value problems, Comput. Math. Appl. 13 (1987) 367–371. J.D. Lambert, Numerical Methods for Ordinary Differential Systems, John Wiley, England, 1991. S.O. Fatunla, Numerical Methods for IVPs in ODEs, Academic Press Inc., London, 1988. S.O. Fatunla, Nonlinear multistep methods for initial value problems, Comput. Math. Appl. 8 (3) (1982) 231–239. S.O. Fatunla, Numerical treatment of singular initial value problems, Comput. Math. Appl. 12B (56) (1986) 1109–1115. F.D. van Niekerk, Rational one-step methods for initial value problems, Comput. Math. Appl. 16 (12) (1987) 1035–1039. S. Abelman, D. Eyre, A numerical study of multistep methods based on continued fractions, Comput. Math. Appl. 20 (8) (1990) 51– 60. M.R. Odekunle, N.D. Oye, s.O. Adee, R.A. Ademiluyi, A class of inverse Runge–Kutta schemes for the numerical integration of singular problems, Appl. Math. Comput. 158 (2004) 149–158. E. Hairer, G. Wanner, Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems, Springer, Berlin, 1996. R.R. Ahmad, N. Yaacob, A.H. Mohd Murid, Explicit methods in solving stiff ordinary differential equations, Int. J. Comput. Math. 81 (2004) 1407–1415. Xin-Yuan Wu, Jian-Lin Xia, Two low accuracy methods for stiff systems, Appl. Math. Comput. 123 (2001) 141–153. J.I. Ramos, Linearization techniques for singularly-perturbed initial-value problems of ordinary differential equations, Appl. Math. Comput. 163 (2005) 1143–1163.
A non-standard explicit integration scheme for initial-value problems Higinio Ramos Scientific Computing Group, Universidad de Salamanca, Spain Escuela Polite´cnica Superior de Zamora, Campus Viriato, 49022 Zamora, Spain
Abstract In this paper we present the construction of a non-standard explicit algorithm for initial-value problems. The method results to be of second order and A-stable. This new algorithm has been proven to be suitable for solving different kind of initial-value problems, specifically, non-singular problems, singular problems, stiff problems and singularly perturbed problems. Some numerical experiments are considered in order to check the behaviour of the method when applied to a variety of initial-value problems. Ó 2006 Elsevier Inc. All rights reserved. Keywords: Non-standard algorithm; Stiff initial-value problems; Singular initial-value problems; Singularly perturbed initial-value problems
1. Introduction A non-standard numerical integration method for initial-value problems (IVP) is presented. It is shown that the new method has second-order accuracy and is A-stable. This favorable property of A-stability suggests that the new scheme may be used for solving stiff systems. Moreover, the non-linear characteristic of the method suggests that it may also be suitable to solve special problems for which conventional numerical integrators result inefficient, such as singular IVPs or singularly-perturbed IVPs. A singular IVP is characterized by the presence of a pole in the solution or by a discontinuous low order derivative of the solution. A simple example is the innocent looking IVP y0 ¼ 1 þ y2;
yð0Þ ¼ 1;
x 2 ½0; b;
ð1Þ
with b > p=4. The theoretical solution is yðxÞ ¼ tanðx þ p=4Þ which has a simple pole at x ¼ p=4. Singularly-perturbed IVPs usually contain a small parameter that multiplies the first-order derivative. These problems are characterized by the presence of thin layers where the solution varies very fast, whereas away from the layers the solution behaves regularly and varies slowly. For the analytical and numerical treatment of this E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.11.134
H. Ramos / Applied Mathematics and Computation 189 (2007) 710–718
711
type of problems the reader is referred to the books by O’Malley [1] and Roos et al. [2]. A common strategy for dealing with this type of problems consists of dividing the domain of integration in two subdomains and then to apply a different scheme on each subdomain (see [3,4]). The use of classical explicit numerical methods for these type of problems requires very small time steps. And the use of standard implicit methods may require a large number of iterations for convergence within each time step. The method presented in this article has been tested on a variety of IVPs and the results show that it may be appropriate for solving IVPs of different nature. The paper has been organized as follows. In Section 2 the non-standard algorithm is developed, in Sections 3 and 4 the main characteristics of the method are established. A section of numerical results on different problems is presented to show the performance of the new method. And, finally, a summary of the main conclusions puts an end to the paper. 2. Non-standard explicit method Consider the IVP y 0 ¼ f ðx; yÞ; yðaÞ ¼ y 0 ; y; f ðx; yÞ 2 R; x 2 ½a; b R;
ð2Þ
where f is assumed to satisfy all the requirements in order that (2) has a unique solution. The interval ½a; b is divided into a number of subintervals ½xj ; xjþ1 with x0 ¼ a and xj ¼ x0 þ jh, being h the step size. Suppose we have solved numerically the problem in (2) up to a point xn, and have obtained a value yn as an approximation of y(xn). Assuming the localization hypothesis [6], yðxn Þ ¼ y n , we are interested in obtaining an approximate value, yn+1, for the true value yðxnþ1 Þ. For that purpose, following the ideas in [5], we suggest an approximation to the theoretical solution yðxn þ hÞ of (2) given by yðxn þ hÞ ¼ yðxn Þ þ
hy 0 ðxn Þ ; aðhÞ þ yðxn Þ
ð3Þ
where a(h) is a sufficiently differentiable unknown function of the step size that has to be determined and it is assumed that aðhÞ þ yðxn Þ 6¼ 0. From (3) it results that F n ðy; a; hÞ ¼ ðyðxn þ hÞ yðxn ÞÞðaðhÞ þ yðxn ÞÞ hy 0 ðxn Þ ¼ 0: Now, if we consider a(h) expanded in Taylor series about h ¼ 0, after expanding F n ðy; a; hÞ in Taylor series about xn the following expression is obtained: 1 F n ðy; a; hÞ ¼ ½1 þ að0Þ þ yðxn Þy 0 ðxn Þh þ ½2a0 ð0Þy 0 ðxn Þ þ ðað0Þ þ yðxn ÞÞy 00 ðxn Þh2 þ Oðh3 Þ: 2
ð4Þ
Imposing that the coefficients of h and h2 in (4) vanish, we obtain a system of equations from which it is readily deduced that að0Þ ¼ 1 yðxn Þ; a0 ð0Þ ¼
y 00 ðxn Þ ; 2y 0 ðxn Þ
provided that y 0 ðxn Þ 6¼ 0. Introducing the above values in the Taylor series of a(h) we have aðhÞ ¼ 1 yðxn Þ
y 00 ðxn Þ h þ Oðh2 Þ: 2y 0 ðxn Þ
ð5Þ
From (3) and (5) it is readily deduced the numerical scheme, which may be written in the form y nþ1 ¼ y n þ
2hf 2n ; 2f n hf 0n
ð6Þ
712
H. Ramos / Applied Mathematics and Computation 189 (2007) 710–718
where y n ¼ yðxn Þ; y nþ1 ’ yðxnþ1 Þ; fn ¼ f ðxn ; y n Þ and fn0 ¼
of of ðxn ; y n Þ þ ðxn ; y n Þfn : ox oy
In particular, when the differential equation in (2) is autonomous, y 0 ¼ f ðyÞ, the method in (6) reads y nþ1 ¼ y n þ
2hfn ; 2 hjn
ð7Þ
ðy n Þ. where we use now the notation fn ¼ f ðy n Þ and jn ¼ of oy Thus we have developed an explicit non-linear one-step method of the form y nþ1 ¼ y n þ hUf ðy n ; xn ; hÞ, where Uf ðy n ; xn ; hÞ is the incremental function, and the subscript f on the right hand side indicates that the dependence of U on yn is through the function f (see [6]). After applying this scheme we will take as approximation for the true solution of (2) at xn+1 the value yn+1 given respectively by (7) or either (6) depending on whether the problem is autonomous or not. Repeating the procedure along the nodes on the integration interval we will obtain a discrete solution for the problem in (2). Non-linear multistep methods are usually designed for dealing with unconventional problems such as stiff or singular initial-value problems, for which the classical schemes generally perform poorly [7,9–12]. On the contrary, some of the non-linear schemes specifically designed for singular problems do not perform well on non-singular problems [12]. In the numerical section we will see that the proposed non-linear scheme is suitable for a variety of IVPs of different nature. 3. Local truncation error It follows from the construction of the method in (6) that the new scheme is at least of second order. In fact, if we consider the functional given by LðzðxÞ; hÞ ¼ zðx þ hÞ zðxÞ
2hz0 ðxÞ2 ; 2z0 ðxÞ hz00 ðxÞ
ð8Þ
where z(x) is an arbitrary function defined on ½a; b and differentiable as often as we need, after expanding in Taylor series about x and collecting terms in h we obtain ! 2 z00 ðxÞ zð3Þ ðxÞ 3 LðzðxÞ; hÞ ¼ þ ð9Þ h þ Oðh4 Þ; 6 4z0 ðxÞ which means that the method has at least second order of accuracy. The local truncation error of the method may be written as ! 2 ðy 00n Þ y 000 n T nþ1 ¼ þ h3 þ Oðh4 Þ; 4y 0n 6 where y 0n ; y 00n ; y 000 n denote, respectively the numerical approximations to the first, second and third derivatives of y(x) at the point xn. Note that if we solve the resulting equation from equating to zero the coefficient of h3 in (9) it may be verified easily that the method in (6) is exact when the solution of the differential equation in (2) is of the form yðxÞ ¼ a þ
b ; cþx
where a; b; c are real constants. This occurs for example for each of the following differential equations: 2
ðaÞ
2y 0 ðxÞy ð3Þ ðxÞ ¼ 3y 00 ðxÞ ;
ðbÞ
2y 0 ðxÞ þ ðc þ xÞy 00 ðxÞ ¼ 0;
ðcÞ
4y 0 ðxÞ þ by 00 ðxÞ ¼ 0;
3
2
H. Ramos / Applied Mathematics and Computation 189 (2007) 710–718
713
2
ðdÞ 2y 0 ðxÞ þ ay 00 ðxÞ yðxÞy 00 ðxÞ ¼ 0; ðeÞ b þ ðc2 þ 2cx þ x2 Þy 0 ðxÞ ¼ 0; ðfÞ yðxÞ þ ðc þ xÞy 0 ðxÞ ¼ a; ðgÞ a2 2ayðxÞ þ yðxÞ2 þ by 0 ðxÞ ¼ 0:
ð10Þ
Remark 1. We observe that when the solution of the differential equation in (2) verifies the equation 2
2y 0 ðxÞ ¼ yðxÞy 00 ðxÞ; then we obtain as a particular case of our method the scheme given by y nþ1 ¼
y 2n ; y n hf ðy n Þ
that is, the first-order scheme in [8], which results to be exact when the solution of the differential equation is of a the form yðxÞ ¼ xþb ; a; b 2 R. This occurs for example in the scheme in (7) for the autonomous problem in the case that y n jn ¼ 2f ðy n Þ holds. 4. Stability analysis The stability of the scheme in (6) is examined by applying it to the Dahlquist’s test problem y 0 ¼ ky; ReðkÞ < 0, yielding the difference equation y nþ1 ¼
2 þ kh y : 2 kh n
Setting z ¼ kh the stability function is obtained as RðzÞ ¼
2þz ; 2z
which is the (1, 1)-Pade´ approximation to the exponential ez, and thus the method is A-stable (see [13]). The stability region of the method consists of the left-half complex plane. 5. Numerical results To test the proposed method, we have integrated several initial value-problems of different nature. All the numerical experiments were performed in double precision using a FORTRAN compiler on a personal computer with a Pentium based processor (2019 MHz). 5.1. A Riccati equation The first example we have considered consists of the Riccati equation coupled with the initial value given by 2
y 0 ðxÞ ¼ 4 þ 4yðxÞ yðxÞ ;
yð0Þ ¼ 1;
whose exact solution is yðxÞ ¼
2x 1 ; x1
and has a pole in x ¼ 1. Although this is a severe difficult problem for most of the usual integrators our method is capable of solving the problem exactly, that is to say, without local truncation error. We have integrated this problem with our method on the interval [0, 2] taking different values for the constant stepsize.
714
H. Ramos / Applied Mathematics and Computation 189 (2007) 710–718
In Table 1 the results for this problem are shown. The errors have been defined as the maximum of the absolute errors on the nodal points in the integration interval Emax ¼ max fjyðxj Þ y j jg; xj 2½a;b
or as the absolute error at the final point Eðx ¼ bÞ ¼ jyðbÞ y NI j; where NI refers to the number of steps used in the integration. The CPU column refers to the computation time measured in milliseconds. We observe that the method has no difficulties in crossing the point x ¼ 1. A plot of the absolute errors is shown in Fig. 1. The differential equation in this problem coincides with the type (g) in (10), and so the errors are only due to the effect of rounding errors. We observe that the bigger the number of steps is so it is the corresponding error. Remark 2. It is possible that for certain combinations of the values appearing in the denominators in (6) or either in (7), some of these denominators vanish. In that case, it is necessary to modify slightly the step size on this step and take mh instead of h, where the factor m is taken nearly unity, say m ¼ 0:95. 5.2. A stiff equation Stiff differential equations involve rapidly increasing or decaying transient solution, which results in a difficulty for most of the numerical integrators. The method in (6) was tested on the stiff problem taken from [14] y 0 ðxÞ ¼ 100yðxÞ þ 99e2x ;
yð0Þ ¼ 0;
which has exact solution yðxÞ ¼
33 2x ðe e100x Þ: 34
The problem has been integrated on the interval [0, 0.5] and the results are presented in Table 2 for different number of steps, NI. Table 1 Errors for y 0 ðxÞ ¼ 4 þ 4yðxÞ yðxÞ2 ; yð0Þ ¼ 1 NI 16 32 64 128 256 512
Eðx ¼ 2Þ
Emax 1:42 10 5:18 1013 2:70 1013 1:99 1011 1:55 1010 1:65 1011
0.00 8:88 1016 3:10 1015 3:55 1015 5:32 1015 2:79 1014
2.5·10
-13
2·10
-13
1.5·10
-13
1·10
-13
5·10
-14
CPU 14
0.5
1
1.5
Fig. 1. Absolute errors for Problem 5.1.
000 000 000 000 000 000
2
H. Ramos / Applied Mathematics and Computation 189 (2007) 710–718
715
Table 2 Errors for y 0 ðxÞ ¼ 100yðxÞ þ 99e2x ; yð0Þ ¼ 0 NI
Eðx ¼ 0:5Þ
Emax
CPU
16 32 64 128 256 512
3:12 102 2:06 105 1:74 106 3:28 107 7:30 108 1:73 108
2:31 101 7:81 102 1:78 102 4:14 103 1:03 103 2:55 104
000 000 000 000 000 016
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5 0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
Fig. 2. Exact and numerical solutions (NI ¼ 64) for Problem 5.2.
In Fig. 2 the exact and the numerical solutions are shown. For NI ¼ 64 we do not appreciate notably differences on the two plots. 5.3. A stiff system The above method may be also applied to a system of equations. If we have y; f ðx; yÞ 2 Rm in (2) we have just to consider the formula in (6) for every component. Let be the stiff system taken from [15] y 01 ðxÞ ¼ 1002y 1 ðxÞ þ 1000y 2 ðxÞ2 ; y 02 ðxÞ
¼ y 1 ðxÞ y 2 ðxÞð1 þ y 2 ðxÞÞ;
y 1 ð0Þ ¼ 1; y 2 ð0Þ ¼ 1;
which will be integrated with our method on the interval [0, 1]. The exact solution is y 1 ðxÞ ¼ e2x ;
y 2 ðxÞ ¼ ex :
The results for every component of the solution appear in Table 3.
Table 3 Errors for Problem 5.3 NI
Eðx ¼ 1Þ y 1 ðxÞ
y 2 ðxÞ
y 1 ðxÞ
y 2 ðxÞ
128 256 512 1024 2048
2:03 103 3:33 105 9:77 108 2:17 108 5:41 109
3:83 103 3:96 104 1:17 107 2:93 108 7:33 109
1:29 102 1:10 102 1:54 107 2:97 108 7:40 109
4:08 103 4:18 104 1:17 107 2:93 108 7:33 109
Emax
CPU
000 000 000 000 000
716
H. Ramos / Applied Mathematics and Computation 189 (2007) 710–718
5.4. A singular problem The next example we consider is the IVP given in (1). The solution becomes unbounded in the neighborhood of the singularity at x ¼ p=4 ’ 0:785398163397448. Finite difference methods based on local polynomial interpolation behave poorly if the solution has singularities. This problem has been used as a numerical test for different non-linear integrators intended for numerically solving singular initial-value problems [7–12]. In Table 4 we present the results for our method where we may observe the ability of the method to cross the singularity. In Fig. 3 the numerical solution after joining the points ðxn ; y n Þ is shown for NI ¼ 128. We observe that the numerical solution reproduces the behaviour of the true solution faithfully. 5.5. A singularly-perturbed problem Now we consider a singularly perturbed IVP given by 2
y 0 ðxÞ ¼ 2kyðxÞ ;
yð0Þ ¼ 10;
which has the exact solution yðxÞ ¼
10 : 1 þ 20kx
We have considered k ¼ 500 as in [16] so that the solution drops quickly from its initial value of 10 to very small values. This problem exhibits an initial layer near x ¼ 0. Table 5 shows the errors obtained with our non-linear method when the integration is performed on the interval [0, 1]. We note that even near x ¼ 0, the new method performs very well, even with the presence of the initial layer there. In fact, the differential equation in this problem coincides with the type (g) in (10), and thus the errors of the method are due only to round off considerations. In Fig. 4 the absolute errors are shown in logarithmic scale when NI ¼ 16.
Table 4 Errors for y 0 ðxÞ ¼ 1 þ yðxÞ2 ; yð0Þ ¼ 1 NI
Eðx ¼ 1Þ
Emax
CPU
32 64 128 256 512 1024 2048
7:18 103 1:79 103 4:48 104 1:12 104 2:80 105 7:00 106 1:75 106
13.9181 3.6385 1.2008 67.1306 16.9897 4.26056 1.06597
000 000 000 000 000 000 016
20 10
0.2
0.4
0.6
0.8
-10 -20 Fig. 3. Numerical solution (NI ¼ 128) for Problem 5.4.
1
H. Ramos / Applied Mathematics and Computation 189 (2007) 710–718
717
Table 5 Errors for y 0 ðxÞ ¼ 2kyðxÞ2 ; yð0Þ ¼ 10; k ¼ 500 NI
Eðx ¼ 1Þ
Emax
CPU
16 32 64 128
2:16 1019 6:50 1019 2:16 1019 0.00
9:71 1017 8:67 1016 7:63 1016 7:21 1016
000 000 000 016
Log (Err) -16 -16.5 -17 -17.5 -18 -18.5 0.2
0.4
0.6
0.8
1
x
Fig. 4. Absolute errors (NI ¼ 16) for Problem 5.5.
5.6. An autonomous problem The last problem to be considered is 3
y 0 ðxÞ ¼ 1 þ yðxÞ ;
yð0Þ ¼ 0;
which will be integrated using different stepsizes in the interval [0, 1.2]. The analytical solution is given implicitly by pffiffiffi pffiffiffi 2y 1 pffiffiffi þ 6 logð1 þ yÞ ¼ 3ð6x þ logð1 y þ y 2 ÞÞ; 3p þ 6 3 arctan 3 from which it is obtained using the Newton method that for x ¼ 1:2 it is yð1:2Þ ¼ 7:368587110472472161317. This value will be used as a reference for calculating the errors. In Table 6 the results obtained with the new method for this problem are presented. In Fig. 5 the numerical solution after joining the points ðxn ; y n Þ is shown for NI ¼ 256.
Table 6 Errors for y 0 ðxÞ ¼ 1 þ yðxÞ3 ; yð0Þ ¼ 0 NI
Eðx ¼ 1:2Þ
CPU
256 512 1024 2048 4096
5:8098 102 1:5670 102 4:0356 103 1:0107 103 2:4024 104
000 000 000 000 016
718
H. Ramos / Applied Mathematics and Computation 189 (2007) 710–718
7 6 5 4 3 2 1 0.2
0.4
0.6
0.8
1
1.2
Fig. 5. Numerical solution (NI ¼ 256) for Problem 5.6.
6. Conclusions In this article the development of a new non-linear explicit scheme for initial value problems is considered. The method is shown to have second order and to be A-stable. These characteristics suggested their use for solving singular and stiff initial-value problems. Numerical results demonstrate that the method can handle these kind of problems at low cost as it is shown by the CPU values in the tables. Acknowledgements The author wish to thank JCYL project SA 024/04 and CYT project BMF 2004-295 for financial support. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]
R.E. O’Malley, Singular Perturbation Methods for Ordinary Differential Equations, Springer-Verlag, New York, 1991. H.G. Roos, M. Stynes, L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations, Springer, Berlin, 1996. S.M. Roberts, A boundary-value technique for singular perturbation problems, J. Math. Anal. Appl. 79 (1982) 489–508. S. Natesan, J. Vigo-Aguiar, N. Ramanujam, A numerical algorithm for singular perturbation problems exhibiting weak boundary layers, Comput. Math. Appl. 45 (2003) 469–479. F.D. van Niekerk, Nonlinear one-step methods for initial value problems, Comput. Math. Appl. 13 (1987) 367–371. J.D. Lambert, Numerical Methods for Ordinary Differential Systems, John Wiley, England, 1991. S.O. Fatunla, Numerical Methods for IVPs in ODEs, Academic Press Inc., London, 1988. S.O. Fatunla, Nonlinear multistep methods for initial value problems, Comput. Math. Appl. 8 (3) (1982) 231–239. S.O. Fatunla, Numerical treatment of singular initial value problems, Comput. Math. Appl. 12B (56) (1986) 1109–1115. F.D. van Niekerk, Rational one-step methods for initial value problems, Comput. Math. Appl. 16 (12) (1987) 1035–1039. S. Abelman, D. Eyre, A numerical study of multistep methods based on continued fractions, Comput. Math. Appl. 20 (8) (1990) 51– 60. M.R. Odekunle, N.D. Oye, s.O. Adee, R.A. Ademiluyi, A class of inverse Runge–Kutta schemes for the numerical integration of singular problems, Appl. Math. Comput. 158 (2004) 149–158. E. Hairer, G. Wanner, Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems, Springer, Berlin, 1996. R.R. Ahmad, N. Yaacob, A.H. Mohd Murid, Explicit methods in solving stiff ordinary differential equations, Int. J. Comput. Math. 81 (2004) 1407–1415. Xin-Yuan Wu, Jian-Lin Xia, Two low accuracy methods for stiff systems, Appl. Math. Comput. 123 (2001) 141–153. J.I. Ramos, Linearization techniques for singularly-perturbed initial-value problems of ordinary differential equations, Appl. Math. Comput. 163 (2005) 1143–1163.