A simple extrapolation method

A simple extrapolation method

21 April 1995 ELSEVIER CHEMICAL PHYSICS LETTERS Chemical Physics Letters 236 (1995)451-456 A simple extrapolation method B,L, Burrows a, M. Cohen ...

512KB Sizes 1 Downloads 262 Views

21 April 1995

ELSEVIER

CHEMICAL PHYSICS LETTERS

Chemical Physics Letters 236 (1995)451-456

A simple extrapolation method B,L, Burrows a, M. Cohen b a Mathematics Group, School of Campmir:g, ,~t,~ffardshireUniversity, Beaconside, Stafford STI80DG, UK b Departmem of Physical Chemistry, The Hebrew University, Jerusalem 91904, Israel

Received 19 October 1994; in final form 24 February 1995

Abstract

We derive and illustrate an extrapolation formula, suitable for monotonic data, which appears competitive with Aitken's classic 8 2 method. A brief error analysis is given which exemplifies the type of sequence for which the formula is most useful. The examples given include sequences of lower bounds to the eigenvalues of one-dimensional SchrSdinger equations. 1. Introduction

There are many problems whose solutions can be obtained only by computing successive approximations and then using some extrapolation algorithm to complete the calculations. Traditional examples include the sums of slowly converging infinite series for which extrapolation formulae use accumulated data on partial sums and iterative solutions of equations of various kinds, but there are many more. In particular, we have been concerned with calculations of sequences of lower bounds to eigenvalues of SchrSdinger's equation [1] and seek to obtain reliable results without resort to laborious calculations. • As is well known, any extrapolation process is generally less reliable than the corresponding interpolation formula, since the extent of the possible error is largely uncontrolled. Nevertheless, there are many algorithms, each based on some assumed plausible functional behaviour, which are employed routinely. Bender and Orzag [2] have reviewed several methods for the treatment of series, while Brezinski [3] has described a general extrapolation algorithm, which includes many of the sequence transforms known to accelerate convergence.

In this Letter, we advocate a particularly simple generalization of the classic 8 2 method due to Aitken (see for example Noble [4]). This formula is usually derived for a sequence {Si} by assuming that Si=S+kt

(1)

i=1,2,3,...,

where S, k and t are constants whose values are found from three successive values S m, Sin+ 1 and Sm+ 2' For convergence to S we require I t ] < 1 and it is clear from (1) that successive errors are assumed to follow the pattern Si+ 1 - S

Si+ 2 - S i + 1 - t =

S, - S

.

Si + 1 - S~

(2)

Despite this restrictive assumption, Aitken's method is used widely and successfully. However there is scope for generalizations of (1) for example by fitting three successive values to alternative forms such as S + a/i Si = 1 + b/--------~'

(3)

where S, a and b are the constants to be found from the data. In Section 2 a general class of extrapolation formulae is considered, which includes the two ex-

0009-2614/95/$09.50 © 1995 ElsevierScience B.V. All rights reserved SSDI 0009-2614(95)00255-3

i,

452

B.L. Burrows, M. Cohen / Chemical Physics Letters 236 (1995) 451-456

amples above and in Section 3 some examples are given. Any discussion of error for extrapolation implies that some assumptions need to be made about the behaviour of the sequence. A brief elementary error analysis is given in Section 4, together with a discussion of the type of sequence appropriate for the extrapolation formula (3).

We assume that we have calculated three approximations S m, S,, Sp which are supposed to be approaching a limiting value S. In most examples m, n and p are consecutive integers, as in the case of successive iterations or in the case of successive partial sums. However, it is perfectly permissible to select any three pieces of data. We adopt the non-linear representation

i - m , n, p,

(4)

v, here S, a and b are constants to be determined by fitting (4) for a suitably chosen function f ( i ) defined on the positive integers. For S i to approach S, f(i) must satisfy the limiting condition

f(i)-.',O

as

i~oo

(5)

and each choice of f(i) is expected to yield a different result in general; these may differ significantly in some cases. Eq. (4) may be rewritten as a set of linear equations for (S, a, b),

S + af( i) - bSif( i) - S,,

(6)

and the solution may be written in the form

SffiSp + a p ( f ) ,

(7)

where

A,(f)

(s,-s.)(s,-s.) = X ( m , n, p ) ( S . - S . , ) - ( S p - S . )

Sm

m
(8)

(10)

(for example, an increasing sequence of lower bounds to an eigenvalue calculated with a variational function that has an increasing number of parameters). Then we expect Ap(f) >i 0 so the choice of the function f(i) must be restricted so that

(S~-S.) X(m, n, p) > (S _Sin ) ,

2. General extrapolation formulae

S + af( i) Sift 1 + b f ( i ) '

For definiteness, suppose that

(11)

which also ensures that the system defined in (6) is non-singular. One possible choice of f(i) is

f ( i ) =t i

(12)

for some real t where ] t ] < 1 so that the limiting condition (5) is satisfied. If we specialise to consecutive integers (m, m + 1, m + 2) and choose

(s.÷2t =

>0,

(Sm+l_Sm)

(13)

then the extrapolation method is identical to Aitken's method defined in (1) with a ffi k and b = 0. In this case 1 X ( m , m + 1, m + 2) = t

(14)

and (7) may be written

S=SA=Sm÷2+

(Sm+2-Sm+l) 2

( s i n ÷ , - s i n ) --(S.,÷2--Sm÷,) '

(15) which is the usual form of Aitken's 8 2 formula, so that the formalism introduced includes Aitken's method. An alternative choice is 1 f(i) - i + c

(16)

where c is any real constant not equal to a negative integer. This leads to p-n

X ( I , n, p ) - - n - m

and

1/f(n)- l/f(p)

X(m, n, p)ffi 1/f(m)- 1/f(n) "

(17)

independent of the choice of c. In particular (9)

X ( m , m + 1, m + 2 ) = 1,

(18)

B.L. Burrows, M. C o h e n / C h e m i c a l Physics Letters 236 (1995) 451-456

which yields

Table 1 Partial sums and extrapolated values for o.j and o"2 (Sm+2-Sm)(Sm+2-Sm+l)

m

S "-- S n - - S m +2 "~" (Sin+l--Sin)

-(Sra+2-Sm+l)

"

(19) This differs from Aitken's formula only in a single factor in the 'correction' term, replacing Sm+ 2 Sm + 1 by the quantity S m + 2 - - Sm" AS we shall see in some examples below, S a is unsuitable for data sequences which are not monotonic. A generalization of (16) is the choice 1

f ( i ) = (i + c) ~'

tr I

0"2

Sm

1 2 3 4 5 6 7 8 9 10

(20)

for a different from unity. Unfortunately, this leads to an expression for X which is not invariant to the choice of c, so that both c and a must be prescribed. However, for m sufficiently large

SA

0.5 0.666667 0.75 0.8 0.833333 0.857143 0.875 0.888889 0.9 0.909091 1.0

m+

2) ~ 1 + ( m + c )

1, m +

n '

(21)

so that X approaches unity asymptotically as In increases. On the other hand if we wish to extrapolate from low values of m, appropriate values of c and a may enhance the rate of convergence. However, the equations for c and a are non-linear if they are determined from additional data.

SB

0.833333 0.875 0.9 0.916667 0.928571 0.9375 0.944444 0.95

Sm

1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0

1.0 1.25 1.361111 1.423611 1.463611 1.491389 1.511797 1.527422 1.539768 1.549768 1.644934

SA

Sa

1.45 1.503968 1.534722 1.554520 1.568312 1.578464 1.586246 1.592399

1.65 1.646825 1.645833 1.645429 1.645235 1.645130 1.645069 1.645031

Table 1 contains a comparison of our results (S a) together with Aitken's (S A) for the series 0-1 and o'2.0-t is clearly a special case, since the partial sum

(a-l) X(m,

453

~n=

1

~ m(m+l) m=l

1

1

= 1 - t~/ ÷ = 1

1 + 1/n

(25)

3. Some examples

is seen to have exactly the assumed form (4); the Riemann zeta function (o"2) is seen to behave quite similarly, and S a is clearly more accurate than S A here. Our comparison for the alternating series for In 2 (o"3) presented in Table 2 suggests that '~'B should not generally be used directly from data which is not monotonic. But for examples such as this, it is a

Slowly converging sequences and series provide a traditional testing ground for extrapolation schemes, and we begin with a few representative examples:

Table 2 Partial sums and extrapolated values for o"3

oo 1 0-1= __~ m ( m + l

)

#t ~-I

1

"*1'0'

'a"2

0"2= )". ~ - - * - - ~ m=l oo

0"3 " -

1.644934,

(23)

6 (_l)m

E

m=l

(22)

-I -'-}

In

In 2 ~ 0.693147.

(24)

Since our aim is to extrapolate from the lowest possible values of m, we confine our attention to a few leading digits only.

m

Sm

1 2 3 4 5 6 7 8 9 10 oQ

1.0 0.5 0.833333 0.583333 0.783333 0.616667 0.759524 0.634524 0.745635 0.645635 0.693147

SA

Sa ~

0.7 0.690476 0.694444 0.692424 0.693590 0.692857 0.693347 0.693003

0.9 0.547619 0.805556 0.601515 0.770513 0.626190 0.752171 0.640372

a Using all the data entries, c Using the even subsequence.

Sa b

SB c

0.690476 0.694444 0.692424 0.693590 0.692857 0.693347

b Using the odd subsequence.

B.L. Burrows, M. Cohen/Chemical Physics Letters 236 (1995) 451-456

454

simple matter to treat one of the monotonic subsequences, with good results as before. A different type of example is the finite-difference solution of an initial value problem, y'=l-x+4y,

y(0)=l.

(26)

Boyce arid Diprima [5, p. 299] present the results of calculations based on Euler's method, using integraI and Y6h. ! tion steps of h ffi 0.1, ~h, ~h To use the . first triple, we assume that mh = ½nh- ¼ph giving X(m, n, p)ffi 2; similarly for the second triple ½mh ffi ¼nh = ~ph giving X(m, n, p ) = 3. Table 3 contains our extra~qlated results, which are seen to provide excellent approximations to the exact solution over the interval 0 < x < 1. However, we note that although the approximate input data consists of a monotonic increasing sequence the extrapolated values overshoot for x > 0.4. In spite of these results, it should not be assumed that S u is universally preferable to SA. When the probable form of the error can be predicted, another algorithm may well be more suitable. For example, if we seek the real root of the cubic equation x 3-3x+

1=0

Table 3 Numerical solution of the initial value problem y' = 1 - x + 4 y , y(O) = 1 x

Extrapolated

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Exact

X=2

X=3

1.60894 2.50507 3.82976 5.79405 8.71307 13.0575 19.5304 29.1828 43.5854 65.0873

1.60903 2.50530 3.83010 5.79421 8.71215 13.0532 19.5174 29.1496 43.5087 64.9210

1.60904 2.50533 3.83014 5.79423 8.71200 13.0525 19.5155 29.1449 43.4979 64.8978

(27)

rationalization of this particular choice is the following. The wave-function parameters in Ref. [1] were obtained by solving a set of m coupled non-linear equations of one general type together with two additional equations, so that c = 2 merely reflects that the true number of parameters is m + 2. Any a greater than unity might be chosen to represent the non-linearity of the problem, and we have made an arbitrary (but apparently highly :.~:ccessful) choice.

(28)

4. Error analysis and discussion

and use iteration based on the scheme xn+ t = ~(1 + x3),

we obtain the results of Table 4. In this case, Aitken's method is predictably preferable, since the errors (8,, 8,+ I) in (x,, x,,+ t) are related by 6,,+ I ~

x28.

In order to obtain a quantitative estimate of the rate o f c o n v e r g e n c e w e c o n s t r u c t a theoretical analytic f u n c t i o n g ( x ) o f the real v a r i a b l e x s u c h that

(29)

and t - x 2 in Eq. (12). Finally, we wish to extrapolate some recently calculated lower bounds to eigenvalues of one-dimensional Schr6dinger equations [1]. The potentials treated were I

V, = - ~.~2 + ~ x

4

,

|

v,. = ~ x

4

,

(30)

appropriate to a symmetric double well and a quartic oscillator, respectively. Our successive approximations are characterized by increasing numbers of variational parameters, and we give eigenvalues for the two lowest-lying even parity solutions in each case; Table 5 contains our results, and in this case only we have designated E c based on the generalization (20) with both c and a chosen equal to 2. Our

g ( x , , ) f f i S n,

1 x,,=-,

n -- I , 2, ... N, ..., M

n

(31)

Table 4 lterative solution of x 3 - 3 x + 1 = 0 R

Xn

1 2 3 4 5 6 7 8 oo

0.4 0.354667 0.348204 0.347406 0.347310 0.347298 0.347297 0.347296 0.347296

XA

X B

0.347130 0.347294 0.347296

0.339593 0.346383 0.347186 0.347283 0.347295 0.347296

B.L. Burrows, M. Cohen / Chemical Physics Letters 236 (1995) 451-456

and so that

g(0)=

lim S,,= T.

the parameters (S, a, b) are chosen by fitting at x N, xN+ 1 and X,v+2 so that

(32)

lq -.i, o o

F,v(xu)=FN(xN+l)=Fu(xs+2)=O

We note that one possible g ( x ) is the interpolation polynomial through the points (1, S 1), (2, $2), . . . ( l / M , SM), (0, T) but any analytical function through these points will suffice. The estimate of error depends on some higher derivatives of g ( x ) so that the optimal choice is an analytic function whose higher derivatives are small. Noise will always be present in the calculated values of S~ from numerical roundoff but we assume that such errors are included in the definition of g(x). However if excessive noise is present then it may be preferable to use an extrapolation procedure that involves least-squares estimates rather than those defined here. We have considered such a scheme elsewhere [6]. For any three successive data points N, N + 1, N + 2 we define the analytic function FN(x)=[l+bf(x)lg(x)-[S+af(x)],

455

and construct the function

GN( X) -- FN( x) -

x 1

1-

XN+ 2

-. XN+ 1

FN(O),

(3Z)

with the properties aN(x,

) = 6 N ( X N + , ) = 6,,(XN+ :) = G,, (0) = 0. (36)

Applying Rolle's theorem to Gu(x), we have G~)(0)=0,

where0<0
u,

(37)

so that, from (35) and (37) e = FN(0 ) = ---~XNXN+IXbI+2F(N3)( , O).

(33)

Table 5 Calculated and extrapolated lower bounds to eigenvalues of Schr6dinger's equation

VI ground state

second exerted state

V2 ground state

second excited state

E.

0 1 2 3 oo a

-0.5 - 0.437994 - 0.432177 - 0.430268 - 0.427588

1 2 3 a

0.472140 0.480051 0.482241 0.485359

a

0.529317 0.529825 0.529966 0.530131 0.530181

a

3.724983 3.726167 3.726679 3.727849

0 1 2 3 1 2 3

(38)

The value of Fu~3)(0) depends on the chosen f ( x ) and the functional form of g(x), which is uniquely defined only at the data points, so that the procedure will yield good results provided that there is some g ( x ) which ensures that Fu~3)(0) is small. The above

where f ( x ) is the extension of f(i) to the real domain (0, ~). We note that FN(0)= g ( 0 ) - S = T - S = e which is a measure of the error when we use the estimate S from the data at N, N + 1 and N + 2. We suppose that, as described in Section 2,

m

(34)

EA

Eu



-0.431502 -0.429419

-0.425009 -0.426721

-0.426714 -0.427590

0.486108

0.483956

0.485199

0.530020 0.530087

0.530215 0.530270

0.530145 0.530156

3.727069

3.727971

3.727609

Ec calculated from (20) with c = ot = 2 (see text). a Variational upper bound accurate to at least the number of digits quoted.

456

B.L. Burrows, M. Cohen / Chemical Physics Letters 236 (1995) 451-456

analysis can be generalized straightforwardly in the case where one (or both) of the non-linear parameters (c, a ) are calculated by fitting at XN+ 3 and/or XN+ 4 leading to an error proportional to N -5 with an unknown coefficient proportional to the fifth derivative of FN(x). To complete the error analysis we need to estimate the derivatives. We note that it i~ also possible to choose c and a in (20) semi-empirically so that F~3)(0) is small and that this is the procedure used for the calculation of E c in Table 5. This procedure, whilst clearly successful in the values obtaiaed in Table 5, is necessarily problem dependent. In the same way certain data will be such that (2) holds approximately and for such problems Aitken's method will be nearly exact and the results will almost certainly be superior to those obtained from (16). An example of such data is given in Table 4. However we can describe a pattern of data for which the latter scheme is well suited. Since c is arbitrary in (16) we choose it to be zero so that f ( x ) - x . For this case

F~( x) - (1 + bx) g( x) - ( S + ax),

(39)

so that

F(~3)(x)-g(3)(x)(1 + bx) + 3bg"(x),

(40)

where b depends on N. For small values of F~3)(0) we require g"(x) and g(3)(x) to be small in the interval 0 < x < 1/N which will be true if we can choose g(x) to be almost linear. In order for g(x) to be approximately linear we require

Xn+ 1 -- Xn

(41)

This is equivalent to n

t=

n+2

< 1

We are grateful to a referee for his valuable comments on an earlier version of this work. The work was initiated during a visit to Oxford by one of us (MC), who enjoyed the generous hospitality of the Department of Theoretical Chemistry and of the Principal and Fellows of St. Edmund Hall, Oxford.

References

(1993) i.

9

.ffi v, N + t . . . . .

Acknowledgement

[I] B.L. Burrows, M. Cohen and T. Feldmann, J. Math. Phys. 34

g(x.÷,)-g(x.) Xa+2 -- Xn+ !

when we have continued, monotonic slow convergence where extrapolation is essential. The example sequences given in Table 1 are typical cases which explains why the results for S B are superior to those obtained using SA. This also explains why in Table 2 the results using either the odd or the even terms are superior to those using all the data, since the former are monotonic sequences. Similarly we can see that the data sets from the variational calculations in Table 5 are monotonically slowly convergent sequences and that the results from S B are better than those using SA. This leads us to believe that the extrapolation procedure defined by (16)will be useful for many sets of data from problems of chemical and physical interest.

(42)

for all n >_.N and t defined by (13). Thus we have good convergence when t = 1 and t < 1; that is

[2] C.M. Bender and S.A. Orzag, Advanced mathematical methods for scientists and engineers (McGraw Hill, New York, 1978). [3] C. Brezinski, Numer. Math. 35 (1980) 175. [4] B. Noble, Numerical methods: iteration, programming and algebraic equations (Oliver and Boyd, London, 1964). [5] W.E. Boyce and R.C. DiPrima, Elementary differential equations (Wiley, New York, 1969). [6] B.L. Burrows and M. Cohen, Chem. Phys. Letters 199 (1992) 580.