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).