European Society of Computational Methods in Sciences and Engineering (ESCMSE) Journal of Numerical Analysis, Industrial and Applied Mathematics (JNAIAM) vol. 5, no. 1-2, 2010, pp. 3-16 ISSN 1790–8140 High Order Finite Difference Schemes for the Solution of Second Order Initial Value Problems 1 Pierluigi Amodio2 , Giuseppina Settanni3 Dipartimento di Matematica, Università di Bari, I-70125 Bari, Italy Received 22 February, 2010; accepted in revised form 13 April, 2010 Abstract: The numerical solution of second order ordinary differential equations with initial conditions is here approached by approximating each derivative by means of a set of finite difference schemes of high order. The stability properties of the obtained methods are discussed. Some numerical tests, reported to emphasize pros and cons of the approach, motivate possible choices on the use of these formulae. c 2010 European Society of Computational Methods in Sciences and Engineering Keywords: Second order Initial Value Problems, finite difference schemes, Boundary Value Methods, absolute stability. Mathematics Subject Classification: 65L05, 65L12, 65L20 PACS: 02.60.Lj, 02.70.Bf 1 Introduction The study of second order differential equations F(t, y, y ′ , y ′′ ) = 0 has a huge bibliography covering several applicative fields, from chemistry to physics and engineering. Even if any high order ODE may be recast as a first order one, this transformation increases the size of the original problem and should make its numerical solution more complicated since it requires the computation of both solution and derivatives (which have different slopes) at the same time. The most interesting and studied second order problems are two-point boundary value problems (see [7] and the reference therein). Among these, two-point singular perturbation problems (see [11]) ǫy ′′ = f (t, y, y ′ ), 0 < ǫ ≪ 1, 1 Published electronically October 15, 2010 author. E-mail: amodio@dm.uniba.it 3 E-mail: settanni@dm.uniba.it 2 Corresponding P. Amodio, G. Settanni 4 have a great appeal since they are stiff, and hence several numerical techniques have been considered for their solution. Among these we recall, for example, collocation methods (see [6]), largely used in the codes. Most of the initial value problems arise from celestial mechanics and lack of the first derivative term. For such problems (called conservative) ad hoc methods have been developed that preserve some properties of the solution [1, 14]. In this paper we are rather interested in initial value problems with nonnull derivative terms since they are not integrated with classical linear multistep formulae. A wide class of such problems arises from singular problems (not defined at some points of the domain, see, for example, [13]) that however will not be considered in this paper. Moreover, some of these are involved in the solution of the nonlinear Schrödinger equation which forms the envelope equation for many physical processes and also for the transverse modulation of a water wave [9, 10]. The idea carried on through this paper is largely used to solve partial differential equations on regular domains. In fact, when it is possible to subdivide the domain with regular grids, each derivative term can be separately approximated. The main gap of this approach is the order of the obtained approximation which is at most 2. In [4] it is suggested how to overcome this problem for BVPs. As a matter of fact, by using the typical approach of BVMs [8] (initial, main and final formulae), it is possible to obtain stable formulae of arbitrary high order. In [5] a generalization of the first order upwind method has been derived for scalar singular perturbation problems. A code based on these formulae is proposed in [2, 3]. Here we apply the same idea to general second order initial value problems. The paper is organized as follows: in the next section we introduce high order finite differences to approximate each derivative of the second order problem. Section 3 concerns with the additional schemes that must be considered in order to use the known value of the first derivative at the initial point. Finally, the last section is devoted to various test examples that are solved by means of both constant and variable meshes. 2 High order finite difference approximations Let us analyze the following second-order initial value problem: ( f (t, y, y ′ , y ′′ ) = 0, y(t0 ) = y0 , y ′ (t0 ) = y0′ , (1) where y0 and y0′ are known values. Let us assume that f is a Lipschitz continuous function in order there exists a unique solution y(t) of the above problem. Let t0 < t1 < · · · < tn (2) be a discretizazion covering all the time-interval or a part of it. The idea proposed in [2, 4, 5] for T the solution of BVPs is that of computing the numerical solution Y = (y0 , y1 , . . . , yn ) of (1) at the meshpoints (2) by approximating y ′ (ti ) and y ′′ (ti ), for i = 1, . . . , n−1, by means of appropriate finite difference schemes (ν) y (ν) (ti ) ≈ yi = k−s 1 X (ν) α yi+j , hνi j=−s j+s (ν) ν = 1, 2, (3) (ν) where hi = ti − ti−1 and the coefficients α0 , . . . , αk are computed such that the formulae have maximum order. In (3) the index k depends on ν and the order of the formula (an order p c 2010 European Society of Computational Methods in Sciences and Engineering (ESCMSE) Finite difference schemes for the solution of second order IVPs 5 approximation for y (ν) (ti ) is in general defined on p + ν points), 0 ≤ s ≤ k and the coefficients depend on k, s, and ν. For Boundary Value Problems, this kind of approach uses variable meshes to discretize the whole time interval. Here, it is preferable to discretize recursively small parts of it by using (for simplicity) constant stepsize inside each subinterval. For this reason, the coefficients of the methods (computed by solving Vandermonde linear systems with integer coefficients) are exactly derived. With respect to the formulae suggested in [2, 4, 5] we have to note that the initial condition y ′ (t0 ) in (1) is not used in (3) which only works with the values yi ≈ y(ti ). Possible approaches to make use of y ′ (t0 ) will be considered in the next section. For the moment, we recall that the methods approximating the derivatives are based on the idea of Boundary Value Methods [8]. For each derivative we fix the order and derive the set of finite difference schemes (3) by changing conveniently the number s and k − s of initial and final conditions, respectively. Among these formulae, we emphasize the main scheme which will be used when possible on the points of the mesh (2). The other formulae (or some of them) will be used once in the extreme points of the mesh. For example, to approximate y ′ (tn ) and y ′ (t1 ) it is necessary to use a final method with only initial conditions and an initial method with at most 1 initial condition, respectively. In vector form the overall approximations for the ν-th derivative can be cast as Y (ν) = 1 Aν · Y, hν ν = 1, 2, where Aν is a (n − 1) × (n + 1) coefficient matrix, whereas Y , containing the n − 1 unknowns of Y , is the solution of the nonlinear system of equations f (Y , Y (1) , Y (2) ) = 0. (4) In the linear case, (4) is a linear system M ·Y =b (5) with the coefficient matrix M having essentially a band structure. We examine only main schemes with approximatively the same number of initial and final conditions. According to several previous papers (see, for example, [4] and the references therein), we call extended central (EC) differences those having the same number of initial and final conditions. The coefficients of EC schemes are symmetric for the second derivative and skew-symmetric for the first derivative. On the contrary, if s again denotes the number of initial conditions, we call generalized backward (GB) and generalized forward (GF) differences those having s − 1 or s − 2 and s + 1 or s + 2 final conditions, respectively, depending on the order of the method. In the following, a subscript after the acronym suggests the derivative to which the scheme is applied. For the second derivative, we consider three possible choices depending on the overall order p. For even orders, we use a symmetric scheme EC2 (defined for k = p and s = k2 ) while for odd orders we choose between generalized forward GF2 or backward differences GB2 (defined for k = p + 1 and s = k2 − 1 or s = k2 + 1, respectively). For the first derivative, we have k = p and s may be k+1 k k k chosen between k−1 2 and 2 if k is odd, and among 2 − 1, 2 , 2 + 1 if k is even. These formulae will be called GF1 , EC1 and GB1 according to what said previously. The combination of such formulae gives rise to 7 couple of main schemes (it is not possible to define symmetric schemes of odd order): EC2 EC1 , EC2 GF1 and EC2 GB1 of even order which were already used for BVPs (see [2, 4, 5]) where it is important to consider a symmetric approximation for the second derivative, and GB2 GB1 , GB2 GF1 , GF2 GB1 and GF2 GF1 of odd order. In the section devoted to the numerical tests, we only consider couple of methods with the same approximation for the derivatives, namely EC2 EC1 , GB2 GB1 , and GF2 GF1 . c 2010 European Society of Computational Methods in Sciences and Engineering (ESCMSE) P. Amodio, G. Settanni 6 As an example, the following are main schemes of order 5 and 6 for the approximation of y ′ and y ′′ . The coefficients of the GB2 (GB1 ) schemes are symmetric (skew-symmetric) with respect to those of the GF2 (GF1 ) schemes. Order 5 GF2 : GF1 : 19 7 10 1 1 1 13 yi−2 + yi−1 − yi + yi+1 + yi+2 − yi+3 + yi+4 180 15 3 9 12 15 90 1 1 1 1 1 ′ h y (ti ) ≈ yi−2 − yi−1 − yi + yi+1 − yi+2 + yi+3 20 2 3 4 30 h2 y ′′ (ti ) ≈ − Order 6 EC2 : EC1 : GF1 : 1 3 3 49 3 3 1 yi−3 − yi−2 + yi−1 − yi + yi+1 − yi+2 + yi+3 90 20 2 18 2 20 90 1 3 3 3 3 1 h y ′ (ti ) ≈ − yi−3 + yi−2 − yi−1 + yi+1 − yi+2 + yi+3 60 20 4 4 20 60 1 2 7 4 1 2 1 h y ′ (ti ) ≈ yi−2 − yi−1 − yi + yi+1 − yi+2 + yi+3 − yi+4 30 5 12 3 2 15 60 h2 y ′′ (ti ) ≈ To investigate the stability properties of the main schemes for the two derivatives, let us analyze their behavior on scalar linear problems of the form y ′′ + γy ′ + µy = 0, (6) where γ and µ are real numbers independent of t. If associated to initial conditions, (6) is well conditioned only when γ and µ are non negative values. In particular, the general solution of (6), y(t) = c1 er1 t + c2 er2 t , where c1 and c2 depend on the initial conditions and r1 and r2 are roots (supposed distinct) of the equation r2 + γr + µ = 0, is monotone decreasing when γ > 0 and bounded for γ = 0. Depending on the initial conditions, y(t) may be strictly positive on all the time interval. Supposing for the moment that the size n of the grid is large, then the effect of the additional methods on the solution may be considered negligible. Therefore, it is sufficient to study the roots of the characteristic polynomial π(z) = ρ(z) + hγσ(z) + h2 µz s̄ , s̄ = max(s1 , s2 ) where ρ and σ are the polynomials associated to the main schemes discretizing, respectively, the second and the first derivative term in (6), and sν , ν = 1, 2, is the value of s in (3) associated to the ν-th derivative. We require that the number of upper off-diagonals of the coefficient matrix in (5) matches the number of roots of π outside the open unit disk [8]. Since γ and µ are real numbers, essentially this means that in the quarter of the plane with hγ ≥ 0 and h2 µ ≥ 0 we have to draw the boundary locus defined by the the straight lines π(1) = 0 and π(−1) = 0, and by the curve π(z) = 0 with |z| = 1 and Im(z) 6= 0. Since z = 1 is always a root of both ρ and σ, the first straight line coincides with the abscissa h2 µ = 0. The condition π(−1) = 0 corresponds to the straight line σ1 hγ + h2 µ ≤ ρ1 , where σ1 and ρ1 are summarized in Tables 1 and 2 for different methods and orders. Since ρ1 > 0, the straight lines corresponding to GB1 , GF1 and EC1 schemes are decreasing (σ1 > 0), increasing (σ1 < 0) and parallel to the hγ-axis (σ1 = 0), respectively. Hence, the use of GF1 schemes give rise to the largest stability domain. Finally, the curve corresponding to the complex values of z of unitary modulus starts from the origin and it is quite near to the h2 µ-axis (it coincides with the segment 0 ≤ h2 µ ≤ ρ1 in case of EC1 schemes). We observe that any combination of formulae does not give stable methods for every value of h (following [8] we say that there is no As,k−s -stable method). As an example, we plot in Figure c 2010 European Society of Computational Methods in Sciences and Engineering (ESCMSE) Finite difference schemes for the solution of second order IVPs 7 Table 1: Coefficients σ1 and ρ1 of the straight line σ1 hγ + h2 µ = ρ1 corresponding to π(−1) = 0. Even order approximations. order 4 6 8 10 EC2 ρ1 16/3 272/45 2048/315 512/75 −64/35 −512/315 GF1 σ1 −8/3 −32/15 GB1 σ1 8/3 32/15 64/35 512/315 0 0 0 0 EC1 σ1 Table 2: Coefficients σ1 and ρ1 of the straight line σ1 hγ + h2 µ = ρ1 corresponding to π(−1) = 0. Odd order approximations. order 3 5 7 9 GF2 /GB2 ρ1 8/3 208/45 352/63 9278/1575 GF1 σ1 −4/3 −16/15 −32/35 −256/315 GB1 σ1 4/3 16/15 32/35 256/315 1 the stability domains for the GF2 GF1 and GB2 GB1 schemes of order 5 and 7. We note that the higher order methods have a larger stability domain. The case γ = 0 is not efficiently solved by this approach when the time interval is very large. In fact, as already known, second order IVPs y ′′ = f (t, y) are properly solved by symmetric linear multistep methods that do not fall in this class of methods. 3 Additional formulae From the number of roots greater than 1 in modulus inside the boundary locus we obtain that, despite the continuous problem has initial conditions, each main scheme (3) requires s − 1 initial and k − s − 1 final formulae. Hence, for example, the EC schemes must always be joined to the same number of initial and final formulae. This section just concerns with the additional schemes we have to use in order to approximate y ′ and y ′′ at both the extreme points of (2). In general we consider formulae in the family (3) with a reduced number of initial or final conditions. Anyway, for second order initial value problems (1), the function f is not approximated at t0 and, hence, the initial value y ′ (t0 ) is not used. Then, we may follow two strategies: the first one is to consider a formula approximating y ′ (t0 ) as an additional equation, thus obtaining the following system  y0 given,    k  1X αj yj = y0′ , (7) h   j=0   f (ti , yi , yi′ , yi′′ ) = 0, for i = 1, . . . , n − 1, T which provides a unique solution Y = (y0 , y1 , . . . , yn ) of the discrete problem. The second strategy is based on the new initial formulae   k X 1 (ν,i) (ν,i) αj yj−1  , ν = 1, 2, y (ν) (ti ) ≈ ν ᾱ0 hy0′ + h j=1 (8) c 2010 European Society of Computational Methods in Sciences and Engineering (ESCMSE) P. Amodio, G. Settanni 8 GF GF schemes of order 5 and 7 2 12 GB GB schemes of order 5 and 7 1 2 10 2 h µ h µ 9 10 1 2 8 order 7 8 7 order 5 6 6 5 4 4 order 5 2 2 0 order 7 3 hγ 0 1 2 3 4 5 6 7 8 9 1 0 hγ 0 1 2 3 4 5 6 7 8 Figure 1: Stability regions for the GF2 GF1 (left) and GB2 GB1 (right) schemes of order 5 and 7. T that use the prescribed value y0′ . Then Y = (y0′ , y0 , y1 , . . . , yn ) is computed by means of  y0 , y0′ given, f (ti , yi , yi′ , yi′′ ) = 0, for i = 0, . . . , n − 1. (9) Both the approaches are applicable to the first block of approximations, corresponding to the initial subinterval. If tn is not the end-point of the time interval (we need other blocks to cover the time domain), from the second block on we could use the last two points as known initial values T of the new block and formulae (3) to uniquely compute Y = (y−1 , y0 , y1 , . . . , yn ) from  y0 , y−1 given, (10) f (ti , yi , yi′ , yi′′ ) = 0, for i = 0, . . . , n − 1. In case the stepsize is changed, the initial value y−1 is computed by means of interpolation techniques from the points in the previous block. As an alternative, it is possible to define a formula analogous to that in (7) to compute an approximation to y ′ (tn ) and then use the same set of formulae considered for the first block. Concerning the approach in (9), for symmetry reason it is more convenient to set Y = (y0′ , y0 , T y1 , . . . , yn , yn′ ) and use (9) (for i = 0, . . . , n) with schemes analogous to (8),   k ν X (−1)  (ν,i) ′ (ν,i) y (ν) (tn−i ) ≈ −ᾱ0 hyn + ν = 1, 2, (11) αj yn−j+1  , hν j=1 as final formulae. We observe that the coefficients in (11) for y (ν) (tn−i ) are just the same as those in (8) for y (ν) (ti ). From a numerical point of view, the last formula, even if it is described in compact form and it is applicable to any block of the time interval, contains values of y and y ′ that could be different in magnitude. Anyway, the numerical tests show that this approach gives the most accurate results. The idea of neglecting the value of the derivative after the first block seems to be more natural for this kind of approach but it requires interpolation formulae that could be ill-conditioned if the order is high since the used stepsize inside each block is constant. Possibly, a variable stepsize inside each block could be advantageous, but we shall not consider this issue here. c 2010 European Society of Computational Methods in Sciences and Engineering (ESCMSE) Finite difference schemes for the solution of second order IVPs 9 As an example, the following are initial formulae in (8) of order 5 and 6 (we recall that y ′ (t0 ) = y0′ ). Order 5 137 ′ 12019 20 5 2 hy0 − y0 + 10y1 − 5y2 + y3 − y4 + y5 30 1800 9 8 25 13 ′ 3281 41 11 5 1 1 2 ′′ h y (t1 ) ≈ hy + y0 − y1 + y2 − y3 + y4 − y5 30 0 1800 12 6 18 24 300 1 1 3 1 1 37 h y ′ (t1 ) ≈ − hy0′ − y0 + y1 + y2 − y3 + y4 4 48 6 4 6 48 h2 y ′′ (t0 ) ≈ − Order 6 49 ′ 13489 15 40 15 hy − y0 + 12y1 − y2 + y3 − y4 + 10 0 1800 2 9 8 77 2171 203 43 13 1 hy ′ + y0 − y1 + y2 − y3 + y4 + h2 y ′′ (t1 ) ≈ 180 0 1200 60 24 54 48 1 ′ 197 1 1 1 1 ′ h y (t1 ) ≈ − hy0 − y0 − y1 + y2 − y3 + y4 − y5 5 300 12 3 12 100 h2 y ′′ (t0 ) ≈ − 4 12 1 y5 − y6 25 18 1 1 y5 − y6 300 1080 Numerical examples In this section we compare both the numerical schemes described in Section 2 and the possible choices for the additional formulae described in Section 3 on two linear and one nonlinear initial value problems. In our numerical experiments we have used blocks with p + 4 equidistant points, where the order p ranges from 3 to 10. For all the examples we have first considered a constant stepsize implementation (see Tables 3, 4, and 5) in order to estimate the order of convergence and to compare the methods. Then, we have solved each problem by means of a simple variable stepsize strategy (see Tables 6, 7, 8, and 9) with initial stepsize h0 = 8 · 10−2 and exit tolerance tol = 10−8 . As usual for this kind of solvers (see [2, 3]), the method of order p + 2 allows us to estimate the error for the method of order p and a new steplength for the successive block by means the formula hnew = 0.9  tol err 1/(p+1) hold . However, in the tables we list the actual absolute error. We focus our attention on three methods which differ each other for the considered additional (initial/final) methods. The first two methods use (7) for the first block. Then the first method (Method 1) computes an approximation for y ′ (tn−1 ) (since the error constant for this formula is much lower than the analogous at tn , the approximation yn at tn is discarded) and iterates on the subsequent blocks. The second method (Method 2) uses yn and yn−1 as initial points in order to apply (10) to the subsequent blocks. Since this approach should require a very carefully interpolation technique, we have not considered it with variable stepsize. The third method (Method 3) uses (9), the initial formulae (8) and the final formulae (11). Problem 1. The first linear problem, y ′′ (t) + y ′ (t) = 0, t ∈ [0, 40], c 2010 European Society of Computational Methods in Sciences and Engineering (ESCMSE) P. Amodio, G. Settanni 10 h = 1 · 10−2 , 4000 points h = 2 · 10−2 , 2000 points h = 4 · 10−2 , 1000 points h = 8 · 10−2 , 500 points Table 3: Numerical results for Test problem 1 with y(0) = 1, constant stepsize. Main Scheme GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 Order 3 3 4 5 5 6 7 7 8 9 9 10 3 3 4 5 5 6 7 7 8 9 9 10 3 3 4 5 5 6 7 7 8 9 9 10 3 3 4 5 5 6 7 7 8 9 9 10 Method1 Error t : y(t) < 0 4.65e-05 4.49e-05 3.31e-05 10.32 2.50e-08 6.65e-08 1.44e-07 15.76 1.23e-10 23.92 1.23e-10 23.92 7.05e-10 21.12 1.15e-12 27.76 6.31e-13 28.40 3.57e-12 26.40 1.54e-05 1.51e-05 3.94e-06 12.48 5.36e-09 6.90e-09 4.36e-09 19.28 2.82e-12 3.02e-12 5.24e-12 26.00 1.90e-13 4.30e-13 28.60 2.30e-13 4.31e-06 4.27e-06 4.79e-07 14.56 4.66e-10 5.14e-10 1.34e-10 22.74 3.53e-13 3.20e-13 28.80 3.74e-14 3.27e-12 7.33e-13 27.96 1.27e-12 1.13e-06 1.13e-06 5.90e-08 16.65 3.36e-11 3.46e-11 4.64e-12 26.10 8.91e-13 27.75 1.14e-12 27.51 4.22e-13 28.57 1.13e-11 1.83e-12 7.58e-12 Method2 Error t : y(t) < 0 3.90e-06 12.96 8.40e-06 11.76 1.89e-05 10.88 6.70e-08 16.56 2.35e-08 17.60 9.03e-08 16.24 2.68e-10 22.08 2.53e-10 22.16 4.80e-10 21.52 1.59e-12 27.20 1.37e-12 27.36 2.60e-12 26.72 1.73e-07 8.06e-07 14.04 1.72e-06 13.28 1.57e-09 20.28 3.31e-10 21.96 2.06e-09 20.04 1.40e-12 27.32 1.43e-12 27.28 2.71e-12 26.64 6.79e-14 3.75e-13 1.17e-13 30.20 4.92e-08 8.36e-08 16.30 1.69e-07 15.60 3.90e-11 23.98 3.72e-12 27.56 4.97e-11 23.74 9.57e-14 2.78e-13 1.69e-13 29.42 1.40e-12 3.38e-13 1.26e-12 8.23e-09 9.32e-09 18.50 1.81e-08 17.83 9.73e-13 1.45e-12 1.05e-12 27.58 2.93e-13 28.89 4.40e-13 28.46 3.62e-13 28.65 3.48e-12 1.61e-12 27.64 4.30e-12 Method3 Error t : y(t) < 0 3.40e-06 6.21e-06 12.00 2.88e-07 16.24 9.86e-09 18.48 1.65e-08 3.57e-09 1.88e-11 1.45e-11 24.96 9.04e-12 2.00e-14 32.56 5.28e-14 3.36e-14 5.13e-07 6.89e-07 14.20 1.78e-08 18.96 3.64e-10 21.76 4.69e-10 5.65e-11 1.62e-13 8.20e-14 30.16 6.20e-14 7.33e-15 1.48e-14 5.53e-14 30.56 6.96e-08 8.06e-08 16.34 1.11e-09 21.74 1.21e-11 25.14 1.41e-11 8.78e-13 2.43e-13 1.67e-13 1.40e-13 9.58e-14 29.98 2.34e-14 31.70 1.41e-13 29.60 9.04e-09 9.73e-09 18.45 6.99e-11 24.58 3.17e-13 1.10e-12 3.03e-13 28.83 6.06e-13 7.39e-13 7.77e-13 4.14e-14 33.80 2.39e-13 29.07 9.27e-13 27.71 c 2010 European Society of Computational Methods in Sciences and Engineering (ESCMSE) Finite difference schemes for the solution of second order IVPs 11 h = 1 · 10−2 , 1885 points h = 2 · 10−2 , 942 points h = 4 · 10−2 , 471 points h = 8 · 10−2 , 236 points Table 4: Numerical results for Test problem 2, constant stepsize. Main Scheme GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 Order 3 3 4 5 5 6 7 7 8 9 9 10 3 3 4 5 5 6 7 7 8 9 9 10 3 3 4 5 5 6 7 7 8 9 9 10 3 3 4 5 5 6 7 7 8 9 9 10 Method1 7.63e-03 7.31e-04 2.46e-03 1.34e-04 1.12e-04 7.66e-05 5.08e-06 5.04e-06 5.15e-06 1.57e-05 1.55e-05 1.62e-05 9.63e-04 1.32e-04 4.03e-04 4.86e-06 4.12e-06 1.34e-06 7.10e-08 6.95e-08 5.18e-08 8.22e-10 6.06e-10 6.62e-10 1.12e-04 2.27e-05 5.56e-05 2.50e-07 2.27e-07 2.93e-08 7.51e-10 8.50e-10 3.95e-10 9.16e-10 9.70e-10 3.97e-10 1.17e-05 5.79e-06 7.26e-06 1.49e-08 1.37e-08 5.10e-10 3.98e-10 9.65e-11 3.32e-10 3.68e-09 3.95e-10 2.76e-09 Error Method2 1.73e-03 6.64e-03 3.83e-03 3.05e-04 2.77e-04 3.35e-04 2.50e-05 2.52e-05 2.59e-05 1.17e-06 1.23e-06 9.85e-07 4.06e-04 6.11e-04 1.82e-04 4.11e-06 3.42e-06 4.30e-06 5.31e-08 5.43e-08 5.15e-08 9.95e-11 6.65e-11 4.85e-10 6.39e-05 6.13e-05 3.69e-06 6.00e-08 4.01e-08 5.99e-08 1.91e-10 1.88e-10 1.82e-10 4.12e-10 1.59e-10 1.57e-10 8.81e-06 6.70e-06 7.76e-07 1.85e-09 1.10e-09 1.09e-09 1.28e-10 3.56e-10 3.25e-10 2.43e-09 1.57e-09 7.96e-10 Method3 2.99e-03 3.30e-03 6.58e-05 3.27e-06 1.38e-05 8.27e-06 9.96e-07 9.32e-07 1.19e-06 3.03e-07 3.28e-07 3.28e-07 3.78e-04 4.00e-04 4.42e-06 1.20e-07 3.39e-07 1.10e-07 3.88e-10 5.25e-10 1.67e-10 7.58e-12 2.09e-11 1.10e-11 4.78e-05 4.91e-05 2.81e-07 5.90e-09 9.14e-09 1.76e-09 3.63e-12 6.02e-11 8.24e-11 1.57e-11 6.85e-12 1.34e-10 6.01e-06 6.09e-06 1.70e-08 6.59e-10 2.17e-10 1.67e-10 3.38e-10 3.04e-10 3.30e-10 6.34e-11 1.77e-10 4.53e-10 c 2010 European Society of Computational Methods in Sciences and Engineering (ESCMSE) P. Amodio, G. Settanni 12 h = 1 · 10−2 , 1000 points h = 2 · 10−2 , 500 points h = 4 · 10−2 , 250 points h = 8 · 10−2 , 125 points Table 5: Numerical results for Test problem 3, constant stepsize. Main Scheme GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 Order 3 3 4 5 5 6 7 7 8 9 9 10 3 3 4 5 5 6 7 7 8 9 9 10 3 3 4 5 5 6 7 7 8 9 9 10 3 3 4 5 5 6 7 7 8 9 9 10 Method1 5.20e-05 3.23e-05 1.10e-04 7.00e-06 6.50e-06 9.29e-06 1.16e-06 1.11e-06 1.36e-06 2.39e-07 2.34e-07 2.68e-07 4.54e-06 3.23e-05 1.31e-05 1.92e-07 1.77e-07 3.31e-07 1.31e-08 1.25e-08 1.68e-08 1.10e-09 1.07e-09 1.29e-09 2.75e-06 3.40e-06 1.56e-06 2.54e-09 2.26e-09 1.03e-08 9.15e-11 8.59e-11 1.46e-10 1.08e-11 5.48e-12 8.68e-12 9.03e-07 9.85e-07 1.89e-07 2.11e-10 2.22e-10 3.11e-10 3.93e-12 3.51e-12 5.28e-12 3.60e-11 5.38e-12 1.38e-11 Error Method2 6.97e-05 3.63e-05 9.47e-05 7.29e-06 6.83e-06 8.91e-06 1.17e-06 1.12e-06 1.35e-06 2.40e-07 2.35e-07 2.68e-07 7.11e-06 2.57e-06 9.19e-06 2.27e-07 2.15e-07 2.85e-07 1.36e-08 1.30e-08 1.59e-08 1.12e-09 1.09e-09 1.26e-09 6.92e-07 1.38e-07 8.35e-07 5.54e-09 5.34e-09 7.19e-09 1.03e-10 9.82e-11 1.21e-10 7.93e-12 6.21e-12 5.87e-12 6.97e-08 7.78e-09 7.82e-08 1.24e-10 1.22e-10 1.66e-10 1.67e-12 2.29e-12 2.29e-12 2.16e-11 1.41e-11 1.18e-11 Method3 1.42e-05 6.13e-06 3.71e-06 7.35e-07 6.34e-07 7.17e-07 5.32e-08 5.57e-08 6.95e-08 1.11e-08 1.12e-08 1.24e-08 1.75e-06 1.23e-06 2.51e-07 1.80e-08 1.20e-08 1.58e-08 4.06e-10 4.62e-10 5.54e-10 3.52e-11 3.53e-11 3.96e-11 2.12e-07 1.79e-07 1.61e-08 3.70e-10 1.58e-10 2.70e-10 6.75e-12 6.46e-12 5.62e-12 1.93e-11 1.92e-11 1.79e-11 2.58e-08 2.37e-08 1.01e-09 7.95e-12 5.75e-12 4.63e-12 3.48e-12 4.38e-12 2.72e-12 1.88e-12 2.13e-12 3.63e-12 c 2010 European Society of Computational Methods in Sciences and Engineering (ESCMSE) Finite difference schemes for the solution of second order IVPs 13 Table 6: Numerical results for Test problem 1 with y(0) = 1, variable stepsize. Main Scheme GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 Order 3 3 4 5 5 6 7 7 8 9 9 10 Method1 Error Mesh 5.35e-06 3.92e-06 4.91e-07 1.69e-08 2.94e-08 3.36e-08 1.11e-09 7.49e-10 3.20e-09 3.45e-10 2.74e-10 5.88e-10 751 898 1038 407 407 434 288 288 288 223 223 223 Method3 Error Mesh 3.03e-07 1.88e-07 8.10e-09 1.61e-08 1.49e-08 4.40e-09 4.93e-10 1.54e-09 5.48e-10 5.39e-11 7.83e-11 8.99e-11 593 817 809 291 331 321 217 193 217 183 183 183 Table 7: Numerical results for Test problem 1 with y(0) = 2, variable stepsize. Main Scheme GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 Order 3 3 4 5 5 6 7 7 8 9 9 10 Method1 Error Mesh 9.36e-06 6.89e-06 1.38e-06 1.25e-08 5.35e-08 1.20e-07 1.04e-08 7.13e-09 1.47e-08 1.82e-09 1.38e-09 2.23e-09 240 275 268 137 137 146 112 112 112 106 106 106 Method3 Error Mesh 7.84e-07 5.62e-07 1.75e-08 3.93e-08 4.34e-08 1.64e-08 4.31e-09 6.83e-09 2.74e-09 1.66e-09 3.18e-09 2.32e-09 169 209 177 101 121 111 97 85 97 85 85 85 c 2010 European Society of Computational Methods in Sciences and Engineering (ESCMSE) P. Amodio, G. Settanni 14 Table 8: Numerical results for Test problem 2, variable stepsize. Main Scheme GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 Order 3 3 4 5 5 6 7 7 8 9 9 10 Method1 Error Mesh 9.09e-05 1.70e-04 1.98e-05 1.52e-06 1.29e-06 8.14e-07 1.08e-06 1.04e-06 7.58e-06 1.45e-07 1.56e-07 3.42e-07 1311 1206 1199 542 560 542 343 343 343 275 275 275 Method3 Error Mesh 2.88e-05 2.88e-05 2.94e-06 4.27e-07 2.52e-06 5.56e-07 1.10e-07 9.34e-08 1.80e-07 1.22e-06 2.16e-06 1.09e-06 1113 1113 761 421 351 371 253 265 253 197 197 197 Table 9: Numerical results for Test problem 3, variable stepsize. Method GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 GF2 GF1 GB2 GB1 EC2 EC1 Order 3 3 4 5 5 6 7 7 8 9 9 10 Method1 Error Mesh 2.57e-06 2.68e-06 7.42e-07 2.65e-08 2.00e-08 6.55e-08 1.60e-08 1.43e-08 1.84e-08 1.18e-08 1.10e-08 1.32e-08 233 247 240 119 119 128 90 90 101 80 80 80 Method3 Error Mesh 2.25e-07 3.58e-07 1.72e-08 7.94e-09 5.62e-09 6.15e-09 1.74e-09 2.07e-09 2.13e-09 2.67e-09 2.77e-09 2.72e-09 193 153 185 101 101 101 73 73 85 71 71 71 c 2010 European Society of Computational Methods in Sciences and Engineering (ESCMSE) Finite difference schemes for the solution of second order IVPs 15 has been solved with initial conditions y(0) = 1 and y ′ (0) = −1 or y(0) = 2 and y ′ (0) = −1. The roots of z 2 + z = 0 are −1 and 0 and, therefore, the exact solution is ye (t) = e−t + c2 , where c2 = y(0) − 1. Tables 3, 6 and 7 are devoted to this example. Tables 3 and 6 solve the problem with the first choice of initial conditions while Table 7 makes use of the second choice. Even if the numerical solution is monotone decreasing, it might tend to a negative value when c2 = 0. For this reason, in Table 3 we also indicate when the numerical solution eventually becomes negative. With constant stepsize we have not observed differences between the two problems (for this reason we have discarded the table associated to the second choice of initial conditions). Vice versa, with variable stepsize, the problem with the second choice of initial conditions has required a much lower number of points. Problem 2. The second linear problem, y ′′ (t) − cos t y ′ (t) + sin t y(t) = 0, t ∈ [0, 6π], has initial conditions y(0) = 1 and y ′ (0) = 1. The exact solution, ye (t) = esin t , has an oscillating solution with period 2π. Problem 3. The nonlinear problem (y(t) + 1) y ′′ (t) − 3(y ′ (t))2 = 0, t ∈ [1, 10], 1 has initial conditions y(1) = 0 and y ′ (1) = − . The exact solution is 2 1 ye (t) = √ − 1. t From all these examples we obtain that the higher order methods compute a more accurate solution also when a large stepsize is used. On the contrary, due to the large size of the systems, with the smallest constant stepsize considered it is not always possible to achieve the best accuracy (see the last block of Tables 3, 4 and 5). Method 2 gives always the worst results. Moreover, its solution for the first problem always gives a negative solution when the stepsize h = 8 · 10−2 . Anyway, the main drawback of all these methods seems to be just that they do not preserve the sign of the solution. 5 Conclusions In this paper we propose to solve second order ordinary differential equations with initial conditions approximating each derivative by means of a set of finite difference schemes. We have derived several methods depending on the choice of the main scheme and the additional formulae. Concerning this last aspect, we have obtained the best results by considering the first derivative at the extreme points (the left one is a known value) inside each block of unknowns. Vice versa we have not observed differences among the possible choices of the main scheme. The choice of the Generalized Forward methods for the first derivative gives larger stability domains and could be more convenient for the most difficult problems. c 2010 European Society of Computational Methods in Sciences and Engineering (ESCMSE) P. Amodio, G. Settanni 16 References [1] R.A. Al-Khasawneh, F. Ismail and M. Suleiman, Embedded diagonally implicit Runge-KuttaNystrom 4(3) pair for solving special second-order IVPs, Appl. Math. Comput. 190 (2007), 1803-1814. [2] P. Amodio and G. Settanni, Variable step/order generalized upwind methods for the numerical solution of second order singular perturbation problems, J. Numer. Anal. Ind. Appl. Math. 4 (2009), 65–76. [3] P. Amodio and G. Settanni, A deferred correction approach to the solution of singularly perturbed BVPs by high order upwind methods: implementation details, in: Numerical analysis and applied mathematics - ICNAAM 2009. T.E. Simos, G. Psihoyios and Ch. Tsitouras (eds.), AIP Conf. Proc. 1168, issue 1 (2009), 711–714. [4] P. Amodio and I. Sgura, High-order finite difference schemes for the solution of second-order BVPs, J. Comput. Appl. Math. 176 (2005), 59–76. [5] P. Amodio and I. Sgura, High order generalized upwind schemes and numerical solution of singular perturbation problems, BIT 47 (2007), 241-257. [6] U. Ascher, J. Christiansen and R.D. Russell, Algorithm 569: COLSYS: collocation software for boundary–value ODEs, ACM Trans. Math. Software 7 (1981), 223–229. [7] U.M. Ascher, R.M.M. Mattheij and R.D. Russell, Numerical Solution of Boundary Value Problems for ODEs, Classics in Applied Mathematics 13, SIAM, Philadelphia, 1995. [8] L. Brugnano and D. Trigiante, Solving Differential Problems by Multistep Initial and Boundary Value Methods, Gordon and Breach Science Publishers, Amsterdam, 1998. [9] C.J. Budd, Asymptotics of multibump blow-up sel-similar solutions of the nonlinear Schrödinger equation, SIAM J. Appl. Math. 62 (3) (2001), 801–830. [10] C. J. Budd, S.Chen and R.D. Russel, New self-similar solutions of the nonlinear Schrödinger equation with moving mesh computations, J. Comput. Physics 1523 (1999), 756–789. [11] J. Cash, BVP Software web page, http://www.ma.ic.ac.uk/∼jcash/BVP software/readme.html. [12] J.R. Cash and M.H. Wright, A deferred correction method for nonlinear two-point boundary value problems: implementation and numerical evaluation, SIAM J. Sci. Statist. Comput. 12 (1991), 971-989. [13] M.S.H. Chowdhury and I. Hashim, Solutions of a class of singular second-order IVPs by homotopy-perturbation method, Physics Letters A 365 (2007), 439-447. [14] Q. Li, X. Wu, A two-step explicit P-stable method of high phase-lag order for second order IVPs, Appl. Math. Comput. 151 (2004), 17-26. [15] R. Lynch and J.R. Rice, A high-order difference method for differential equations, Math. Comp. 34 (1980), 333–372. c 2010 European Society of Computational Methods in Sciences and Engineering (ESCMSE)