VSVO multistep formulae adapted to perturbed second-order differential equations

Appl. Math. Lett. Vol. 11, No. 3, pp. 83-87, 1998 Copyright©1998 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0893-9659/98 $19.00 + 0.00

Pergamon

PIh S0893-9659(98)00037-8

V S V O Multistep Formulae A d a p t e d to P e r t u r b e d Second-Order Differential Equations J. VIGO-AGUIAR Department of Pure and Applied Mathematics, University of Salamanca E-37008 Salamanca, Spain

jvigo@gugu, asal. es

J. M. FERR/~NDIZ Department of Applied Mathematics, Universidad de Alicante E-03080 Alicante, Spain (Received April 1997; accepted June 1997)

Communicated by D. G. M. Anderson A b s t r a c t - - W e introduce and study multistep formulae of arbitrary order with step variation facilities. They are adapted to integrate second-order differential equations with constant coefficients, in the sense that the homogeneous IVP is integrated exactly (with only the round-off error), and the perturbed problems with a local error that comes only from the perturbing terms. Keywords--Numerical integration of perturbed ODEs, Oscillatory solutions, Variable step-size algorithms.

1. I N T R O D U C T I O N There exists a variety of methods which have been adapted for numerical solution of the IVP (a,/7 real),

L(x) - x " + # x ' + a x = ¢(x, x', t), x(o)

=

t e [0, b],

(1)

=

that all have the property of integrating the homogeneous equation L ( x ) -- 0 exactly (i.e., with only round-off errors). Usually they are restricted to the case/3 = 0, and have been designed to obtain very accurate solutions in satellite orbital dynamics [1,2], highly oscillatory electric circuits [3], and in the Schr6dinger equation [4]. Most of these methods do not allow step size variation, though there are a few exceptions due to Deuflhard [5], Denk [3], and Simos and Raptis [4]. However, these exceptions actually require a quite difficult computation of coefficients, and no multistep or Runge-Kutta algorithm greater than sixth order has been derived so far. By the time that Stiefel and Bettis set up their famous methods [6], a colleague of theirs, Scheifele, had developed a one-step procedure by approximating the solution to (1) as a truncation of the exact one, in the form n

=

bkC

(t).

(2)

k=O Typeset by ~ 4 ~ T F ~ 83

84

J. VIGO-AGUIAR A N D J.

M. FERRANDIZ

The functions Gk are defined as solutions to the IVPs: tn--2 L ( G n ) = (n - 2)!'

Gn(0) = 0,

G'(0) = 0,

L(Go) = 0,

G0(0) = 1,

G~(0) = -f~,

L(G1) = 0,

GI(0) = 0,

G~(0) = 1,

n > 2,

(3)

and may be calculated efficiently using some recursion formulae [7]. The coefficients bk are

bo = xo,

bl = X~o + ~xo,

bk+2 = f(k)(O),

(4)

with f ( t ) = ¢(x(t), x'(t), t) and x(t) the solution to (1). They are calculated using recurrence relations obtained as in the Taylor integration method. Scheifele recommended changing the order rather than the stepsize, this giving the algorithm great accuracy, although it is impossible to apply unless the function ¢ is simple. In the case ~ = 0, Martfn and Ferr£ndiz overcame this difficulty by transforming that method into a fixed-step multistep one [8], whereas the authors in [9] developed the first high-order multistep method that allows feasible step and order variations, since the coefficients are given by explicit expressions and can be calculated easily. In this paper, we introduce the corresponding VSVO scheme for the general second-order equation (1), with the help of some auxiliary results relative to divided differences that are collected in Section 2.

2. S O M E R E L A T I O N S I N V O L V I N G D I F F E R E N C E S Given an increasing sequence of points ti in [a, b] and choosing n, let us denote

hi = ti - ti-1 > 0,

Hi = ti - tn-1,

(5)

so that Hn = hn , H n - 1 = O, H n - 2 = - h n - 1 , H n - 3 = - ( hn-1 + h n - 2 ) , . . . , and let f [tn , tn-1, . . . , tn-k] be the kth-order divided difference of f ( t ) in the values tn, t n - X , . . . , t,,-k. For the sake of brevity, the difference t J - l [ H n , . . . ,Hn-i+l] of a power function will be denoted by q i j ( n ) . We also introduce the k-vectors

Dk,n = (f[tn], f[tn, t n - 1 ] , . - . , f [ t n , . . . , tn-k+l]), ( ft(tn-1) f(k-1)(tn-1)) Fk,n= f(tn-1), 1! ''"' -~:~)v. '

(6)

and the k x k matrix Pk,n = (qi,j(n)). We have the following theorem. THEOREM 1. If H = m a x { H n , . . . ,

H n - k + l } , then D kt , n = Pk,n " F k,n t + Rtk

with Rk = ( O ( H k ) , O ( H k - 1 ) , . . . ,

(r)

O(H)).

PROOF. If f is holomorphic on [a,b], using the Taylor expansion at tn-1 of f ( t n ) . . . . , f(t,~-i+l) leads us to the expression oo

f[tn, tn--l,..., tn--i+l] = Z Cj f(j) ( t n - 1 ) ~=o

J!

(s)

for certain constants Cj that only depend on the H~'s and not on f. To compute Cj, we make tn-1 = 0 and replace f[tn, t n - 1 , . . . , tn-i+l] by tk[Hn, 0 , . . . , Hn-i+l] in (8). As (tk)(J)(0) = k! 6k,j (the Kroneeker delta), we get Cj = tJ[Hn, O,... ,Hn-i+l] = qi,j+l(n). Truncation of (8) at the

VSVO Multistep Formulae

85

term j = k - 1 produces an error of order k - i + 1, since qi,j+l(n) has order j + 1 - i in H (see [10, p. 7])• If f is not analytic, we can follow a similar procedure using finite Taylor expansions instead of (8). | In [11], it is shown that the relation (7) can be solved explicitly for Fk,n. Let a i j ( n ) be the i.e., the sum of all products of j - i distinct variables H n - i , so that ai,i = 1 and

jth elementary symmetric polynomial in H n , . . . , g n - j + 2 , ai,j(n) = ( - 1 ) j - i

Z

(9)

H n - k , Hn-k2 "" " H n - k j _ ,

O
and Sk,n = (ai,j(n)) a k x k matrix, where we have defined a i j ( n ) to be zero for j < i. According to [11] we have the following theorem.

THEOREM 2. T h e matrices Pk,n = ( q i j ( n ) ) and Sk,n are inverses of each other. Moreover, F k,n t = Sk,n " D k,n t + Qtk

(10)

with Qtk of the same order in H that Rtk.

3. T H E

NUMERICAL

METHODS

Taking into account (4), and the original approximation (2) of Scheifele, there evolves a pair of formulae Xn = X n - l ( G o ( h n ) + ~ G l ( h n ) ) + Xln_:G:(hn) + Wk,n " F~, Xln = X l n _ l G o ( h n )

- Xn-l(aGl(hn))

(11)

+ IfVk,n " F~,

where Wk,n = (G2(hn), 1! G 3 ( h n ) , . . . , (k - 1)[ Gk+l(hn)),

(12)

~Yk,n = (Gl(hn), I[ G 2 ( h n ) , . . . , (k - 1)! Gk(hn)).

T h e y provide a multistep integration procedure suitable for equation (1), that can be expressed in terms of differences according to Theorem 2. DEFINITION 1. T h e S V F m u l t i s t e p schemes are de/~ned to be the following: * explicit scheme: xn = X n - : ( G o ( h n ) + ~ G : ( h n ) ) + x~n_:V:(h~) + Wk,n

•

Sk,n-1

D kt, n - 1 ,

"

(13)

x~ = X~_lGo(h~ ) - X n - l ( a G : ( h n ) ) + Wk,n " Sk,,~-: " D kt , n - 1 , • implicit scheme: t xn = x,~-:(Go(h~) + ~3G:(h~)) + x ~ _ : G : ( h n ) + Wk+l,n " Sk+x,,~ " Dk+:,~,

(14)

x~ = X ~ _ l G o ( h n ) - x , ~ - : ( a G : ( h ~ ) ) + l?dk+l,~ • Sk+l,n" D kt+ l , n "

Both schemes have an implementation similar to that of Adams methods with variable coefficients (see, for instance, [12]). In each step, the coefficients of Sk,,~, Wk,,~, and l)¢'k,,~ can be calculated from recurrences. More precisely, we have the following. PROPOSITION

1.

(a) Since ~i,j(n) = a ~ - l , j - l ( n ) - H n - j + l a i , j - l 0 1 Hn o

(n), the m a t r i x Sk,n can be obtained as follows:

0

...

a2,5(n)

""

0

\ 0

1

,

)

:2,3(~)

a2,4(n)

\ Sk,,~ =

~

\ ~

> a2,k(n)

\ aa,5(n)

"'"

'

a3,k(n)

'

a4,k(n)

0

0

a3,4(n)

0

0

1

a4,5(n)

"'"

:

:

:

:

"..

:

0

0

0

0

.-.

1

\

(b) G,~(h) + ~ G n + l ( h ) + aGn+2(h) = (hn/n[).

~

\

,

(15)

86

J. VIGO-AGUIARAND J. M. FERP~NDIZ

PROOF. (a) is an immediate consequence of the definitions, and (b) is a result of Scheifele [7]. The main properties of SVF algorithms are given in the following theorem. THEOREM 3. If the function ¢(u, v, t) of equation (1) is Lipschitzian in the variab/es u and v, then both methods (13) and (14) are convergent 1 of order 2 k and k + 1, respectively. The principal terms of the local error are

k Gj+l(hn)uj,k+l+a(n - 1) f [ t n - 1 , - . . , t n - k - 1 ] ,

for the function x,

j~-O k

(16)

y ~ Gj(hn)~j,k+l+a(n -- 1) f [ t n - l , . . . , tn-k-ll, j~-o

for the derivative x',

where a = 0 for the explicit scheme and a = 1 for the implicit. PROOF. We use the fact that consistency and stability imply convergence (see [13, Theorem 1.2.4, p. 13]). To obtain consistency it is sufficient to consider (7) and to take into account that Gi(hn) has order i in H . The stability is derived from Grigorieff [14, p. 362, Theorem 1]: since the last terms in equations (13) and (14) satisfy a uniform Lipschitz condition, the problem reduces to the stability of the linear part, with ¢ = 0, that is obvious. REMARK 1. The schemes exactly integrate L(x) = O, since here the local truncation error vanishes. If ¢ is considered as a small perturbation, i.e., it has a small parameter as a factor, the local truncation error contains the same factor. These are the reasons why the previous schemes are quite useful for the accurate integration of (1). REMARK 2. We could have used modified divided differences (see [12, p. 143]) with minimal changes in the formulation. Our methods when k = 1 and k = 2 are not used in simulations because their accuracy is too low for the purposes for which the methods are usually employed. The explicit method for k = 2 is presented here merely to illustrate the structural aspects of the procedure. Let us denote

Gi =- Gi(hn) and ~hn 2 Go --- e 5" cOSWn,

G2 = wn2 + 6-------~ 1 - Go +

G1 = e g" sinwn,

G1

,

(17)

hn

The explicit method (13) can be expressed as

Xn ~-~Xn-l(Go + j3G1) "4-x'n_lG1 -~- G2f(tn-1) + Gzf[tn-l,tn-2], / Xn ~---Xtn_lC0 --Xn--laCl Jr" Clf(tn-1) + C2f[tn-1, tn--2].

(zs)

Examples of the efficacy and efficiency of an early version of this method in some Celestial Mechanics model problems can be found in [9]. In addition, in that paper, the eighth- and tenthorder explicit methods are compared to standard and nonstandard codes like the one proposed by Deuflhard [5], a symplectic integrator and its constant step predecessor SMF. XThe definitionsof consistency, stability,and convergence used appear in [13]. 2The order is k ifthe principalterm of the local truncation error is O ( H k+2) for the function x and O(H k+1 ) for the derivativect.

VSVO Multistep Formulae

87

REFERENCES 1. E. Stiefel and G. Scheifele, Linear and Regular Celestial Mechanics, Springer, Berlin, (1970). 2. J. Vigo-Aguiar, Mathematical methods for precise calculation of satellite orbits, Ph.D. Dissertation, University of Valladolid, (1993); (Spanish language, available from the author). 3. G. Denk, A new numerical method for the integration of highly oscillatory second-order ordinary differential equations, App. Num. Math. 13, 57-67 (1993). 4. T.E. Simos, A new variable-step method for the numerical integration of second-order IVP and their application to one-dimensional SchrSdinger equation, Appl. Math. Lett. 6 (3), 67-73 (1993). 5. P. Deuflhard, A study of extrapolation methods based on multistep schemes without parasitic solutions, ZAMP 30, 177-189 (1979). 6. E. Stiefel and D.G. Bettis, Stabilization of Cowell's method, Numer. Math. 13, 154-175 (1969). 7. G. Scheifele, On the numerical integration of perturbed linear oscillating systems, ZAMP 22, 186-210 (1971). 8. P. Martin and J.M. Ferr~ndiz, Multistep numerical methods based on Scheifele G-functions with application to satellite dynamics, SIAM J. Numer. Anal. 34, 359-375 (1997). 9. J. Vigo-Aguiar and J.M. Ferr~ndiz, High-order variable step algorithms adapted to the integration of perturbed oscillators, Computer Physics Communications (submitted). 10. L.M. Milne-Thomson, The Calculus of Finite Differences, The Macmillan Press, New York, (1981). 11. J. Vigo-Aguiar, An approach to variable coefficient multistep methods for special differential equations, Discussion Paper VA197, Dept. Matem£tica Aplicada, Univ. Alicante, (1996). 12. J.D. Lambert, Numerical Methods for Ordinary Differential Systems, John Wiley and Sons, New York, (1991). 13. H.J. Stetter, Analysis of Discretization Methods for Ordinary Differential Equations, Springer, Berlin, (1973). 14. R.D. Grigorieff, Stability of multistep methods on variable grids, Numer. Math. 42, 359--377 (1983).