Noname manuscript No. (will be inserted by the editor) Two-step fourth order methods for linear ODEs of the second order Mikhail V.Bulatov · Guido Vanden Berghe Received: date / Accepted: date Abstract A family of two step difference schemes of the fourth order has been developed for linear ODEs of the second order. Stability proporties for such schemes are discussed and results of numerical tests are given. Keywords Difference schemes · Numerov method · P -stability 1 Introduction In this work new approaches are considered for the development of numerical methods for the solution of problems of the form y ′′ (t) = A(t)y(t) + f (t), y(0) = y0 , y ′ (0) = y0′ , t ∈ [0, 1] , (1) where A(t) is a (n × n) matrix, f (t) is a n-dimensional vector-function. From here on it is supposed that elements of the matrix A(t) and the vector-function f (t) have enough smoothness. Let’s set an equal mesh on a segment [0, 1] ti = ih, i = 0, 1, . . . , N, h = 1/N and let’s denote with yi an approximation of y(ti ) , Ai = A(ti ), fi = f (ti ). Classical k−step methods for the problem considered have a form k X j=0 αj yi+1−j = h2 k X βj (Ai+1−j yi+1−j + fi+1−j ). j=0 M. V. Bulatov Institute of System Dynamics and Control Theory, Sibirian Branch of the Russian Academy of Science, Irkutsk, Russia E-mail: [email protected] G. Vanden Berghe Department of Applied Mathematics and Computer Science, University of Gent, Gent, Belgium Tel.:+32-9-2644805 Fax: +32-9-2644995 E-mail: [email protected] 2 If αj = αk−j and βj = βk−j , j = 0, 1, . . . , [k/2] such methods are called symmetric methods. It is known (see for example,[1],[2]), that the order of a stable k−step method does not exceed k + 2 if k is an even number and k + 1 if k is an odd number.When k = 2 one of the well known methods is the Numerov method yi+1 − 2yi + yi−1 = h2 /12[Ai+1 yi+1 + fi+1 + 10(Ai yi + fi ) + Ai−1 yi−1 + fi−1 ], which is of fourth order. A vast bibliography regarding the numerical solution of the problems considered is given in the monographs [1], [2] and in the paper [3]. In this paper we shall consider the development of a non-classical one-stage two-step difference schemes for the problem (1). Under certain conditions we shall demonstrate that the constructed schemes have order 4. Similar schemes for ODEs of order 1 have been considered in [4],[5]. 2 Two-step schemes Let us consider the following two-step scheme for the problem (1) yi+1 − 2yi + yi−1 = h2 [β01 (Ai+1 yi+1 + fi+1 ) + β11 (Ai yi + fi ) + β21 (Ai−1 yi−1 + fi−1 )] . (2) This method is stable and it has at least first order if it satisfies the condition (see, for example [1],[2],[6]) (3) β01 + β11 + β21 = 1 . By analogy with scheme (2) we can write its one leg variant, i.e. yi+1 − 2yi + yi−1 = h2 [A(β02 ti+1 + β12 ti + β22 ti−1 )(β02 yi+1 + β12 yi + β22 yi−1 ) + +f (β02 ti+1 + β12 ti + β22 ti−1 )] . (4) Coefficients β02 , β12 and β22 will satify a condition analogous to (3), i.e. β02 + β12 + β22 = 1 . (5) Taking into account condition (5) and the fact that ti+1 = ti + h scheme (4) can be re-written in the form yi+1 − 2yi + yi−1 = h2 [A(ti−1 + (β12 + 2β02 )h)(β02 yi+1 + β12 yi + β22 yi−1 ) + f (ti−1 + (β12 + 2β02 )h)] .(6) In order to simplify the notation we denote ti−1 + (β12 + 2β02 )h = tB , remembering that this value changes from step to step. We re-write the relations (2) and (6) in the form of the overdetermined system of linear algebraic equations: E − h2 β01 Ai+1 E − h2 β02 AB −E + h2 β21 Ai−1 −E + h2 β22 AB yi+1 = 2E + h2 β11 Ai 2E + h2 β12 AB yi−1 + h2 ψi φB , yi + (7) where E is the nth-order identity matrix, ψi = β01 fi+1 + β11 fi + β21 fi−1 and φB = f (tB ) = f (ti−1 + (β12 + 2β02 )h), AB = A(tB ). System (7) does not necessarily have 3 a solution in the classical sense; therefore we construct a “solution” as follows. We multiply both sides by the n × 2n rectangular matrix Li+1 = (E − h2 β01 Ai+1 d(E − h2 β02 AB )), where d is a free parameter. In this way we obtain a two-step difference scheme of the form Mi+1 yi+1 = Pi yi + Qi−1 yi−1 + gi+1 , (8) where Mi+1 = Li+1 Qi−1 = Li+1 E − h2 β01 Ai+1 E − h2 β02 AB ; Pi = Li+1 2E + h2 β11 Ai 2E + h2 β12 AB ; (9) −E + h2 β21 Ai−1 −E + h2 β22 AB ; gi+1 = h2 Li+1 ψi φB Since the schemes (2) and (6) have at least first order when (3) and (5) are satisfied, also scheme (8) shall have at least first order for any values of the parameters {βji }, i = 1, 2; j = 0, 1, 2, satifying (3) and (5) and for any value of d. It means that (8) has five free parameters, which can be determined by order conditions, i.e. we develop (8) in series with respect to orders of the step length h. The first order conditions are given by (3) and (5); by this the method has immediately also order two. The order three condition reads (10) d + 1 = β11 + 2β01 + dB, B = β12 + 2β02 . There are two fourth order conditions: 5 (d + 1) = β11 + dβ12 6 5 (d + 1) = β11 − dB 2 + 2dB 6 (11) (12) The fifth order conditions read: Bd(−2 + 3B − B 2 ) = 0, 5 1 1 + B − B2 + B3) = 0 d(1 + d)( 12 12 2 (13) (14) Let us first consider the solutions of the equations making all coefficients in the series expansion of (8) up to powers h4 zero, i.e. solving the system of equations: β01 + β11 + β21 = 1 β02 + β12 + β22 = 1 d + 1 = β11 + 2β01 + dB, 5 β11 + dβ12 5 6 (d + 1) = 1 (d + 1) = β − dB 2 + 2dB 1 6 This system has the following solutions d = 0, β11 = 5 , 6 β01 = 1 , 12 β21 = 1 , 12 (15) 4 free parameters β12 , β02 , β22 , B and 1 1 1 2 3 B − B + 1, β02 = B 2 − B, β12 = −B 2 + 2B, 2 2 2 2 5 β11 = (1 + d) + dB 2 − 2dB, 6 1 11 1 1 1 1 3 1 β0 = (1 + d) − dB 2 + dB, β21 = − d + − dB 2 + dB, 12 2 2 12 12 2 2 β22 = free parameters B, of the form (16) d 6= 0. Remark that solution (15) results in a difference equation 1 2 1 2 5 h A(t + h))2 y(t + h) − (E − h A(t + h))(2E + h2 A(t))y(t) − 12 12 6 1 2 1 2 h A(t + h))(−E + h A(t − h))y(t − h) − (E − 12 12 1 2 1 5 1 (E − h A(t + h))h2 ( f (t + h) + f (t) + f (t − h)). (17) 12 12 6 12 (E − 1 2 After division by the common factor (E − 12 h A(t + h)) it represents nothing else than the well-known Numerov method, of which the truncation error is O(h6 ). The solution (16) has still two free parameters and gives rise to the following difference equation. 1 1 1 (d + 1) − dB 2 + dB)A(t + h))2 + 12 2 2 1 2 1 2 d(E − h ( B − B)A(t − h + hB))2 )y(t + h) = 2 2 1 1 5 2 1 ((E − h ( (d + 1) − dB 2 + dB)A(t + h))(2E + h2 ( (d + 1) + dB 2 − 2dB)A(t)) 12 2 2 6 1 1 −d(E − h2 ( B 2 − B)A(t − h + hB))(2E + h2 (−B 2 + 2B)A(t − h + hB)))y(t) + 2 2 1 1 ((E − h2 ( (d + 1) + dB(1 − B))A(t + h)) × 12 2 1 3 2 1 (−E + h ( (1 − 11d) − dB 2 + dB)A(t − h)) − 12 2 2 1 2 1 2 d(E − h ( B − B)A(t − h + hB)) × 2 2 3 2 1 2 (−E + h ( B − B + 1)A(t − h + hB)))y(t − h) + 2 2 1 1 1 1 (E − h2 ( d + − dB 2 + dB)A(t + h)) 12 12 2 2 1 1 2 1 1 2 − dB + dB)f (t + h) − ×h (( d + 12 12 2 2 5 11 1 1 3 2 ( (d + 1) + dB − 2dB)f (t) + (− d + − dB 2 + dB)f (t − h)) + 6 12 12 2 2 1 2 2 1 2 dh (E − h ( B − B)A(t − h + hB))f (t − h + hB). (18) 2 2 ((E − h2 ( The truncation error of (18) is O(h5 ). Since we still have two free parameters (i.e. B and d) it is still meaningfull to consider equations (13) and (14). The solution of this extended system is: d = 0, β11 = 5 , 6 β01 = 1 , 12 β21 = 1 , free parameters β02 , 12 β22 , β12 , B, 5 which delivers again the Numerov method, d = −1, β12 = 0, B = 0, β11 = 0, β01 = 0, b20 = 0, β22 = 1, β21 = 1 β02 = 1, β22 = 0, β21 = 0 , and d = −1, B = 2, β12 = 0, β11 = 0, β01 = 1, which both are reproducing Mi+1 = Pi = Qi−1 = gi+1 ≡ 0 and which are of no realistic use and B = 1, β12 = 1, 1 5 1 (d + 1), β11 = − d, 12 6 6 β02 = 0, β22 = 0, β01 = β21 = 1 (d + 1), 12 (19) with free parameter d 6= −1. For this last solution the truncation error is O(h6 ). It is not possible to choose for d an appropriate value, which increases this order. 3 Stability properties The linear stabilty theory for the considered methods is descrived by Lambert and Watson [7] and is related to the test equation y ′′ = −k2 y . For this test equation one obtains, with ν = hk yi+1 − 2Rmm (ν 2 )yi + yi−1 = 0 . where the rational function Rmm (ν 2 ) is called the stability function of the method. The interval of periodicity of the method is (0, ν02 ) if |Rmm (ν 2 )| ≤ 1 for 0 ≤ ν 2 ≤ ν02 . The method is P -stable if |Rmm (ν 2 )| ≤ 1 for all real ν. Let us remark that the Numerov method possesses the highest order of all two-step methods, but is has a narrow stability interval (0, 6). In this section we study the stability properties of method (8) with the coeeficients (19) delivering a method of order four depending on one free parameter d 6= −1. This scheme has the following general structure h2 (d + 1)Ai+1 )2 + dE)yi+1 12 h2 h2 −((E − (d + 1)Ai+1 )(2E + (10 − 2d)Ai ) + dE(2E + h2 Ai ))yi 12 12 h2 h2 (d + 1)Ai+1 )(−E + (d + 1)Ai−1 ) − dE)yi−1 (20) −((E − 12 12 h2 1 = h2 (E − (d + 1)Ai+1 )( ((d + 1)fi+1 + (10 − 2d)fi + (d + 1)fi−1 ) + dfi ) 12 12 ((E − and corresponding stability function R22 = 1− 1+ ν2 3 ν2 6 4 ν − (5 − d) 144 4 ν + (d + 1) 144 . (21) 6 Table 1 The end point of the stability intervals for different values of the parameter d d ν02 0 0.5 1 1.5 2 √6 −4 + 4√7 = 6.58... −6 + 6 √ 5 = 7.41... −12 + 12 3 = 8.78... 12 One can easily check that |R22 (ν 2 )| ≤ 1 iff (4 − 2d) ν2 ν4 + −2≤0. 144 6 For d ≥ 94 this condition is always fulfilled and the corresponding methods are P -stable. For d < 49 the methods have a rather large interval of periodicity. In table 1 we give some examples. Remark that for d = 0 we obtain the same results as for Numerov’s method. Graphics of the functions (21) for different values of d are shown at figures 1 and 2. The pictures confirm the theoretical results discussed above. Fig. 1 Stability functions of formula (21) for various d. 4 Numerical experiments In this section we give numerical results for two test examples obtain by using scheme (20) for different values of parameter d. We assume that when realising the difference schemes we know the exact value of the input data y0 and y1 . 7 Fig. 2 Stability functions of formula (21) for various d. Example 1. Consider the equation y” = a2 2a + 4 3 t t y, Z e−2a/x dx, which has a general solution a/t y(t) = C1 e a/t + C2 e where C1 and C2 are arbitrary constants. Let C1 = 1, C2 = 0, x ∈ [1, 10], a = −20. Introduce on the segment [1,10] a mesh with with a step hj = 0.2/2j , j = 1, 2, 3 and let erj = |e−2 − yNj |, Nj = 9/hj . The results of calculations are given in 4. Table 2 d -2 0 2 3 5 h = 0.1 6 · 10−3 3 · 10−4 5 · 10−3 7 · 10−3 1 · 10−2 h = 0.05 6 · 10−4 3 · 10−5 5 · 10−4 8 · 10−4 1 · 10−3 h = 0.025 6 · 10−5 2 · 10−6 4 · 10−5 6 · 10−5 1 · 10−4 Example 2. Use the following problem as a test: ′′ ′ x = −k2 x + 2ak cos(kt), x(0) = 1, x (0) = 0 ′′ ′ y = −k2 y + 2ak sin(kt), y(0) = 0, y (0) = k − a, t ∈ [0, 40π], which has the exact solution of the form x = cos(kt) + at sin(kt), y = sin(kt) − at cos(kt). 8 When k = 1 and a = 0.0005 we get a known test example from [8]. We shall consider the values k = 5 in our numerical experiment. By analogy with [3] we introduce the function γ(t) = and its difference analogue p x2 (t) + y 2 (t) = γhj (i) = p 1 + (at)2 q x2i + yi2 , where xi , yi are found using schemes (20). The method (20) with the input data defined by makes it possible (for a certain value of the parameter d) to take a considerable bigger integration step. √ Also note that when d = 1 the difference schemes proposed are stable for ν 2 ∈ (0, 6( 5 − 1)), and when d = 2 for ν 2 ∈ (0, 12) and P −stable for d ≥ 49 . Introduce on the segment [0, 40π] an analytical grid with the step hj = π/4j , j = 1, 2, 3, 4 erj = |γ(40π) − γhj (Nj )| where Nj = 40 ∗ 4j ∗ hj . Results of the computations for various d and for the steps hj = π/4j , j = 1, 2, 3, 4 when k = 5 are given in table 3. Table 3 The errors for the numerical experiments for k = 5. The symbol * denotes errors larger than 106 . d 1 2 3 4 5 j=1 ∗ ∗ 2.9 · 10−2 4.0 · 10−2 1.5 · 10−1 j=2 1.8 · 10−2 4.6 · 10−2 1.9 · 10−2 7.9 · 10−4 1.7 · 10−2 j=3 9.8 · 10−5 3.4 · 10−4 6.0 · 10−4 8.0 · 10−4 1.0 · 10−3 j=4 3.7 · 10−7 1.2 · 10−6 2.3 · 10−6 3.2 · 10−6 4.4 · 10−6 The results confirm the stability regions and the P -stability properties as obtained theoretically. 5 Conclusions Thus, the schemes proposed possess a number of advantages comparing to the schemes developed previously. Firstly, similarly to the Numerov methods they have order 4 but their stability intervals are bigger than that of the Numerov method for small values of the occurring parameter d. Once the parameter d exceeds the values 9/4 the methods are even P -stable. Secondly, when applying the Obrechkoff methods it is necessary to compute derivatives of the input data which is an ill-posed problem. The proposed methods do not require the computation of derivatives. However, we have to compute the matrix A2 (t) at each step of integration. In some cases we can find the matrix A2 (t) in advance; this is the case if the matrix coefficients are polynomials or if the matrix A(t) has a special structure. The authors believe that investigating properties of non-symmetrical methods and constructing similar schemes for a more general case of problem y ′′ (t) = A(t)y ′ (t) + B(t)y(t) + f (t), can be promising. y(0) = y0 , y ′ (0) = y0′ , t ∈ [0, 1] 9 Acknowledgements The work has been supported by Russian Foundation for Basic Research (project No. 07-01-90000), Registered project No. 5 of Siberian Branch of Russian Academy of Sciences and Funding For Scientific School No. 1676.2008.1 sponsored by President of Russia. References 1. E. Hairer, S. P. Norsett, G. Wanner, Solving Ordinary Differential Equations 1, Nonstiff Problems, Springer-Verlag, Berlin Heidelberg New York, 528p. (1987). 2. L. Gr. Ixaru , G. Vanden Berghe, Exponential Fitting, Kluwer Academic Publisher, Dordrecht,Boston,London, 308p., (2004). 3. M. Van Daele, G. Vanden Berghe, P-stable Obrechkoff methods of arbitrary order for second-order differential equations, Numer. Algor., 44,115-131, (2007). 115-131. 4. M. V. Bulatov, Construction on a one-stage L-stable second order method, Differential equations, 39, 593-595 (2003). 5. M. V. Bulatov, On construction of nonclassic difference schems for ordinary differential equations. Differential equations, 44, 546-557 (2008). 6. N. S. Bakhvalov, Numerical Methods. Moscow ’Nauka’,(1975) (in russian). 7. J. D. Lambert, I. A. Watson, Symmetric multistep methods for periodic initial value problems, J. Inst. Math. Appl., 18, 189-202 (1976). 8. T. E. Stiefel, D. G. Bettis, Stabilization of Cowell’s method, Numer. Math., 13, 154-175 (1969).

© Copyright 2019