- Email: [email protected]
Chapter IV: Applications
195
.
CBAPTER N:Applications
5 1. Convergence 1.1. 1.2. 1.3. 1.4. 1.5.
The The The The The
acceleration . . . . . . . . . . . . . . . . . . . . univariate t-algorithm . . . . . . . . . . . . . . . qd-algorithm . . . . . . . . . . . . . . . . . . . . algorithm of Bulirsch-Stoer . . . . . . . . . . . . . p-algorithm . . . . . . . . . . . . . . . . . . . . . multivariate €-algorithm . . . . . . . . . . . . . .
197
. 197 205
. 209 213
. 216
$2. Nonlinear equations . . . . . . . . . . . . . . . . . 2.1. Iterative methods based on Pad4 approximation . 2.2. Iterative methods based on rational interpolation 2.3. Iterative methods using continued fractions . . . 2.4. The qd-algorithm . . . . . . . . . . . . . . . 2.5. The generalized qd-algorithm . . . . . . . . .
. . . . . 220 . . . . . . 220 . . . .
. . . . . 227 . . . . . 233 . . . . 233 . . . . . 236
53. Initial value problems . . . . . . . . . . . . 3.1. The use of Pad4 approximants . . . . . 3.2. The use of rational interpolants . . . . 3.3. Predictor-corrector methods . . . . . 3.4. Numerical results . . . . . . . . . . . 3.5. Systems of first order ordinary differential
. . . . .
. 238 . 241 . 243
. . . . .
. . . . .
. . . . .
. . . . .
. . . . . . . . . . . . . . . . . . . . equations . . . .
238
244
. 247
.
54 Numerical integration . . . . . . . . . . . . . . . . . . . . . 250 4.1. Methods using Pad4 approximants . . . . . . . . . . . . . 251 4.2. Methods using rational interpolants . . . . . . . . . . . . 252 4.3. Methods without the evaluation of derivatives . . . . . . . . 253 4.4. Numerical results for singular integrands . . . . . . . . . . 254
55 . Partial differential equations
. . . . . . . . . . . . . . . . . . 257
§6.Integral equations . . . . . . . . . . . . . . . . . . . . . . . 260 0.1. Kernels of finite rank . . . . . . . . . . . . . . . . . . 260 6.2. Completely continuous kernels . . . . . . . . . . . . . . 262
Problems
. . . . . . . . . . . . . . . . . . . . . . . . . . .
265
Remarks
. . . . . . . . . . . . . . . . . . . . . . . . . . .
267
. . . . . . . . . . . . . . . . . . . . . . . . . .
268
References
196
“It i s my hope that by demonstratdng the e m e with which the various transformations may be effected, their field of application might be widened, and deeper insight thereby obtained into the problems for whoae aoh6tion the tran8formatdons have been uaed.
P . WYNN (1956).
- “On a
device for computing the e m ( & )
tran~f~rmation~’
Z V . l . l . The univariate €-algorithm
197
The approximations introduced in the previous chapters will now be used to develop techniques for the solution of various mathematical problems: convergence acceleration, numerical integration, the solution of one or more simultaneous nonlinear equations, the solution of initial value problems, boundary value problems, partial differential equations, integral equations, etc. Since these techniques are based on nonlinear approximations they shall be nonlinear thernselves. We shall discuss advantages and disadvantages in each of the sections separately and illustrate their use by means of numerical examples.
§I. Convergence acceleration. 1.1. The univariate r-algorithm. Consider a sequence { a ; } ; , ~ of real or complex numbers with lim ai = A
i-
00
Since we are interested in the limiting value A of the sequence we shall try to construct a sequence (b;),.,N that converges faster to A , or
We shall describe here some nonlinear techniques that can be used for the construction of { b;},N. Consider the univariate power series
i-t
with Vai = a; - a;-1. Then clearly for the partial sums
we have
Fk(1)= ak
k = 0, 1,2, . . .
If we approximate f(z) by r ; , ; ( z ) , the Pad4 approximant of order (i,i) for f , then we can put b; = r;,;(l)
i = 0 , 1 , 2 , . ..
198
I V . l . l . The univariate €-algorithm
For the computation of b; the c-algorithm can be used:
Then
The convergence of the sequence {biIiEm depends very much on the given sequence { a , } ; , N . Of course the convergence properties of { b ; } ; , ~are th e same as those of the diagonal Pad6 approximants evaluated a t z = 1 and for this we refer to section 4 of chapter 11. In some special cases it is possible to prove acceleration of the convergence of ( a ; } i E N . A sequence {a;}iEN is called totally monotone if
Ak u ; > O
d,k=0,1,2
,...
where A k a , = Ak-' ai+i - Ak-' ai and Ao a; = a;. In other words, {ai}iEN is totally monotone if
a0 2 a1 2 a2 2 ... 2 0 Aao 5 Aal 5 h a 2 5 . .. 5 0 A2ao 2 A 2 a l 2 A2a2 2 . . . 2 0 and so on.
IV.1.1. The uniuariate 6-algorithm
199
Note t ha t every totally monotone sequence actually converges. The sequences (Xi}i,-~ for 0 <_ X 5 1 and {l/(i l ) } ; , ] ~are for instance totally monotone sequences. The close link with the theory of Stieltjes series becomes clear in the following theorem [4 p.811.
+
Theorem 4.1. The sequence { o ~ } ~ €isN totally monotone if and only if there exists a realvalued, bounded and nondecreasing function g on [0, 11 taking on infinitely many different values such that 1
a; = I 0
t' d g ( t )
i = 0,1,2,. ..
A sequence { a ; } ; , ~is called totally oscillating if the sequence
ai}iE]N
is totally monotone. One can see th at every convergent totally oscillating sequence converges to zero. For these sequences the following results can be proved [4 pp. 83-85]. Theorem 4.2.
If the 6-algorithm is applied to the sequence { U ; } ~ ~ ] N with 1imi-m a i = A and if there exist constants a and @ such that { a a ; + @ } i E ~is a totally monotone sequence, then lim
cfi
=A
lim
($2
k >_ 0 and fixed
=A
e 2 0 and k e d
f-00
k+m
If also lirn
i-oo
then
ai+l
-A # 1
-_-
ai-A
200
IV.1.1. The univariote t-algorithm
Theorem 4.3.
If the €-algorithm is applied t o the sequence { a ; } i E mwith limi-a; = A = 0 and if there exist constants cr and B such that (&a; + @ ) i E ~is a totally oscil-
lating sequence, then
We illustrate these theorems with some numerical results. Consider
This is a totally monotone sequence. In table 4.1. we have listed the values for n = 0 , . . ., 4 and m = 0, . . ., 10. So the values b i = can be found on the main diagonal. Clearly the convergence of the sequence { a i } i e ~is accelerated.
EE)
All these computations were performed using double precision arithmetic (56digit binary representation). To illustrate the influence of data perturbations and ( m - 4 obtained using single rounding errors we give in table 4.2. the same values E~~ precision input and single precision arithmetic (24-digit binary representation). In the tables 4.3. and 4.4. the respective results are given for the totally oscillating sequence
Here the convergence of the sequence {b;}+-N is much faster th a n th a t of the sequence {ai}i,-N. These examples illustrate that the influence of data perturbations and rounding errors can be important. When is nearly equal to one can loose a lot of significant digits. In [62] i t is shown how this misfortune may be overcome.
tP1)
Table 4.1.
56-digit binary reprea ent otio n %
m\n
S 0
1
2
3
4
0
0.10000000D+01
0.86666667D+00
0.52173913D+OO
0.43636364D-bOO
0.37874803D-I-00
1
0.50000000D+00
0.25000000D+00
0.17847058D+00
0.13801453DfOO
0.113Q8283D-tOO
2
0.33333333D-tOO
0.16668667D+00
O.llllllllD~00
0.85271318D-01
0.69561659D-01
3
0.25000000D+00
0.12500000D+00
0.83333333D-01
0.62500000D-01
0.50607287D-01
4
0.20000000D~00
0.10000000D+00
0.66666667D-01
0.50000000D-01
0.40000000D-01
5
0.16666667D+00
0.83333333D-01
0.55555556D-01
0.41666667D-01
0.33333333D-01
6
0.14285714D-kOO
0.71428571D-01
0.47610048D-01
0.35714286D-01
0.2857142QD-01
7
0.12500000D+00
0.62500000D-01
0.41666667D-01
0.31250000D-01
0.25000000D-01
8
O.lllllIllD~OO
0.55555556D--01
0.37037037D-01
0.27777778D-01
0.22222222D-01
9
O.lOOOOOOOD+OO
0.50000000D-01
0.33333333D-01
0.25000000D-01
0.20000000D-01
10
0.9090Q091D-01
0.45454545D-01
0.30303030D-01
0.22727273D-01
0.18181818D-01
k !Y
% F:
2. c P
3.
R
cs
hl
0
N
Table 4.2.
24 -digit binary repreaent ation
m\.
0
1
2
3
4
c..
0
0.10000000E+01
0.86666680E+00
0.52173018E+00
0.436363703+00
0.378748OOE+OO
k
1
0.500000003~00
0.250000083+00
0.17647058E+00
O.l3801455E+OO
0.113082853~00
Y
2
0.33333334E-kOO
0.16666664E+00
0.11111123E+00
0.852712323-01
0.69561742E-01
E
3
0.25000000E+00
0.12500004E+00
0.83333202E-01
0.62500395E-01
0.50605543E-01
4
0.20000000E+00
0.10000000E+00
0.66666812E-01
0.49907089E-01
0.40013745E-01
5
0.16666667E+00
0.83333351E-01
0.55554848E-01
0.416750723-01
0.332807823-01
6
0.14285715E+00
0.71428463E-01
0.47620941E:-01
0.356953373-01
0.28692538E-01
7
0.12500000E+00
0.62500104E-01
0.41664604E:-01
0.312740473-01
0.247051133-01 0.22477847E-01
8
O.llllllllE~00
0.55555556E-01
0.370379803-01
0.27756441E-01
0
0.10000000E+00
0.5000001OE-01
0.33332549E-01
0.2502411OE-01
O.l0564494E:-Ol
10
0.90909004E-01
0.45454517E-01
0.30304417E-01
0.226919743-01
0.18706985E-01
z 2. C
e
2.
f: m
m
it
rg 0 7
3
Table 4.3.
56-digit binary repreaentatran Ls,
m\.
5
0
0
0.10000000D+01
1
-0.50000000D~00
2
0.33333333D.fOO
3
4 5 6
7 8
-0.25000000D+00
1 0.40000000D~00 0.35714286D-01 -0.98039216D-02 0.40322581D-02
0.20000000D~00 -0.20408163D-02 -0.16666667D+00 0.14285714D-kOO -0.12500000D+00 0.111111llD+OO
0
-0.100OOOOOD+00
10
0.~0009001D-01
0.1 17370881)-02 -0.73637703D-03 0.40212508D-03 -0.34506556D-03 0.25125628D-03 -0.18860807D-03
2 0.25531015D-kOO -0.2 18078 10D-01 0.1 1454754D-02 -0.25960540D-03 0.83229297D--04 -0.330401 llD-04 0.15170805D-04 -0.77490546D-05 0.42861980D-05 -0.252506 123)-05 0.1 5651583D-05
3 0.18604651 D+OO 0.11833001D-01 -0.38107008D-03 0.34806706D-04 -0.73003358D--05 0.20314308D-05 -0.688Q106OD-06 0.27006940D-06 -0.11826320D-06 0.58510073D-07 -0.289770043)-07
4
0.14600870D+00 -0.96160448D-02 0.15341245D-03 -0.89086733D-05 0.10459978D-05 -0.20044207D-06 0.53279358D-07 -0.16203117D-07
; Y
z C
2. C P -.
l
a
cs m
a
cq
*
a 0
0.5656006OD-08 -0.220288523)-08 0.0378304OD-00
ts 0 W
01
0
Table 4.4.
b b
m+l Rd-digit binary representation
m n 0
0 0.10000000E~01
1
-0.500OOOOOE~00
2
0.333333348+00
3
4 5
-0.25000000E+00 0.20000000E+00 -0.16666667E+00
6
0.142857153+00
7
-0.12500000E~00
8
9 10
O.llllllllE~00 -0.10000000E+00 O.QO909094E-01
1
0.400000043+00 0.357142843-01 -0.880381823-02 0.40322733E-02 -0.20408179E-02 0.11736Q87E-02 -0.738376743-03 0.48213349E-03 -0.345061333-03 0.25125383E-03 -0.18860842E-03
2 0.25531918E-kOO -0.218878OOE-01 0.1145486OE-02 -0.25958874E-03 0.832260823-04 -0.3305387QE-04 0.15181173E-04 -0.77439190E-05 0.428871413-05 -0.25255961E-05 0,156517563-05
4
3 0.18604654E3+00 0.118331053-01 -0.38106500E-03 0.34898013E-04 -0.73038805E-05 0.20300561E-05 -0.68641026E-06 0.273580183-08 -0.11687172E-08 0.56633620E-07 -0.28493082E-07
0.146908823+00 -0.9816957OE-02 0.15340716E-03 -0.881250783-05 0.10442981E-05 -0.20828127E-06 0.56101026E-07 -0.138258973-07 0.860953913-08 -0.17431754E--08 0.183647753-08
* Y
k k
Y
z E
a.
a 3. C
h R
m
$-
-l
3
I V . 1.2. The qd- algorithm
205
Various generalizations t o accelerate the convergence of a sequence of vectors, matrices, and so on exist. We refer to [4] and [Sl]. The convergence of a multidimensional table of real or complex numbers is treated in section 1.5. 1.8. The qd-algorithm.
To accelerate the convergence of { ai}iEN we construct again
i- 1
We have seen in section 3 of chapter 11 th at it is possible t o construct a corresponding continued fraction I
I
I
where the coefficients d i can be calculated by means of the qd-algorithm. If we calculate the i t h convergent Ci(z) of this continued fraction, we can now put 6i = Gi(1)
i = 0 , 1 , 2 , .. .
The numbers 6; can be computed using a continued fraction algorithm from section 5 of chapter I. Remark that theoretically the numbers b 2 ; here are the same as the 6; computed by means of the €-algorithm. In practice there are however differences due to rounding errors. We will illustrate this numerically. The tables 4.5. and 4.6. illustrate the effect of perturbations when the qd-table is computed: the results in table 4.5. are obtained using double precision arithmetic while those in table 4.6. are obtaincd using single precision arithmetic. Input was
So f(z)= ez and limi-m a; = e = 2.718281828.. . From the tables 4.5. and 4.6. the values b; were computed using the three-term recurrence relations (1.3.) respectively in double and single precision. The results can be found in the respective tables 4.7. and 4.8., where they can also be compared with the values of the c-algorithm.
Table 4.5.
hl
0 Q,
k=l,2,
92
-0.16668667D-kOO 0.33333333D+00
-0.23800524D-01 0.142857 14D+00
-0.17857143D-01 0.12500000D+00
0.97222222D-01
-0.1388888OD-01 0.1 11111 11D+OO
-0.90900091D-02 0.90000001D-01
0.727272731)-01
0.18181818D-01
0.81818182D-01
-0.27272727D-01
-0.45454545D-01
0.54545455D-01
-0.36363636D-01
0.63636364D-01
(4
e5
0.55555556D-01
-0.44444444D-01
0.66666667D-01
-0.33333333D-01
0.77777778D-01
-0.22222222D-01
0.8888888OD-01
-0.11111111D-01 0.10000000D+00
-0.27777778D-01
(kf
95
-0.55555556D-01
0.60444444D-01
-0.41666667D-01
0.83333333D-01
(a)
e4
0.71428571D-01
-0.53571429D-01
0.80285714D-01
-0.35714286D-01
0.107 14286D-kOO
(k)
94
-0.7142857lD-01
0.05238005D-01
-0.47610048D-01
0.110047821)+00
(k) e3
0.10000000D+00
-0.66666667D-01
0.13333333D-kOO
-0.33333333D-01 0.16866867D+00
(k)
93
-0.10000000D+00
0.15000000D~00
-0.50000000D-01 0.20000000D+00
(k)
e2
O.I6666667D+OO
-0.83333333D--01 0.25000000D+00
56-digit binary reprerent ation
(k)
0.50000000D~00
...
Y t r
.4
n
Table 4.6. gd - table f o r gik’ =
Vak+1
v ak
_ I _
k
with ak =
C -j!1
k = 1,2, ...
j-0
4 - d i g i t binary representation
0.333333343+00
-0.833333433-01 0.25000000E~00
0.16666666E+00
0.14285713E-l-00
-0.17857134E-01 0.12500000E~00
-0.13888896E-01 0.111 llllOE+OO
0.888887723-01
-0.11111103E--01 0.100000003+00
-0.454895578-01
0.545119843-01
-0.36346123E-01
0.836478143-01
-0.272756813-01
0.72725780L?3-01
-0.1 81818443-01
0.818182303-01
-0.0090007OE-02 0.80809084E-01
-0.222224373-01
(k)
e5
0.555540553-01
-0.444474973-01
0.666824538-01
-0,333314843-01
0.77779025E-01
(4
95
-0.55559866B--01
0.694431893:-01
-0.416667613-01
0.83332881E-01
-0.277776273-01
0.972223 13E-01
(k)
e4
0.714307433-01
-0.535723123-01
0.89285374E-01
-0.357142693-01
0.107142813$00
(i)
94
-0.714278823-01
0.052385073-01
-0.47619178E-01
O.l1804755E+OO
-0.238085223-01
(4
e3
0.99Q89800E-01
-0.666665883-01
0.13333338E+00
-0.333333483-01
93
-0.10000007E+00
O.l4999998E+OO
-0.49999989E-01 0.20000000E+00
(k)
0.16666670E+00
IV.1.2. The qd-algorithm
208
Table 4.7.
lim a; = 2.71828182845904..
i-+co
j-0
56-digit binary repreaentation (0)
i
€2 i
1.000000000000000D+00 1.000000000000000D+OO 3.000000000000000D+OO 3.000000000000000Df00 2.7142857142857 14D+00 2.714285714285714D+OO 2.7 18309859154930D+00 2.718309859 154930D+00 2.718281718281718D+OO 2.718281718281718D+OO 2.7 18281828735696D+00 2.718281828735696D+OO 2.7182818284585631)+00 2.7 18281828458564D+00 2.718281828459046D+OO 2.718281828459046D+OO 2.7 18281828459045D+00 2.718281828459045D+OO Table 4.8. 1
a i = j -C 0,
lim ai = 2.718281828.. .
i-
00
24-digit binary represent ation 1
0 1
2 3 4 5
I 1.0000000E+00 3.0000000E+00 2.7142858E+00 2.7183099E+00 2.7 1828183+00 2.7182815Ei-00
1.0000000E+00 3.0000000E+00 2.7 I42858E+OO 2.71830993+00 2.71828183+00 2.71828223+00
.
IV.1.3. The dgorithm of Bulirach-Stoer
209
1.3. The algorithm of Bulirsch-Stoer. Let
(2,)
i
e be ~ a convergent sequence of points with
limi-+w z i = z When using extrapolation techniques to accelerate the convergence of we compute a sequence {bi};,N with
{U;},.~N
b; = limz-.z gi(z) where g i ( z ) is determined by the interpolation conditions gi(zi) = a j
j
0 , . . ., 9: with i
2
0
The point z is called the extrapolation point. Often polynomials are used for the interpolating functions 8 ; . We consider here the case of rational interpolation where we shall use Stoer’s recursive scheme given in section 3 of chapter 111. Using the notations of chapter I11 and writing
we first generalize the formulas (3.12.):
and
with pO,o(s) ( i ) = a, and qo,o(z) (i) = 1. We can easily see taking into account the interpolation properties of r(J+l)
m,n- 1
rk!n and
that
(4.1a .)
IV.1.3. The algorithm of Bulirech-Stoer
210
and analogously that
(4.1b .)
If aiT;ll
# 0, then we can write
or
Analogously
If we use the interpolating functions
IV.f.3.The algorithm of Bulirach-Stoer
and write 8 (; )' = r@,(z) for rn
with
= t-tb(z)= a j and
211
+ n = &,then
8-1 ( j1
= 0 for j = 0 , 1 , 2 , . . , and
Let ua compare t h e algorithm of Bulirsch-Stoer with the €-algorithm, the qdalgorithm and the method of Neville-extrapolation based on polynomial interpolation (here bi = rri(z)). Input is the sequence {ai);,N = { ~ / r n ) ; ~ N
with lim;+a a; = 0. For the methods based on interpolation we use
z=
lim x i = O
i- oa
IV.i.3. The algorithm of Bulirsch-Stoer
212
The results are listed in table 4.9. Computations were performed in single precision arithmetic.
Table 4.9. a.=
1
-__
&Ti Bulirsch-
i ~
qd
Neville
Stoer
epsilon
Qi _____--_-- b i
s(0)
€ i(0)
C;(1)
0
1.00000000
1.00000000
1.00000000
1.00000000
1.00000000
2
0.5773502e
0.26964909
0.24118084
0.47414386
0.47414389
4
0.4 472 13 59
0.16024503
0.11 8 2 0 6 8 3
0.3095461 1
0.3095461 1
6
0.3779e450
0.11422185
0.07439047
0.22860453
0.22960511
8
0.33333334
3.08880618
0.06251375
0.18284784
0.18285738
10
0.30151 1 3 5
0.07302472
0.05196425
0.15728275
0,15730451
12
0.27735010
0.0617 1018
0.04684320
0.14 1 I5 857
0.14100026
14
0.25819880
0.11071037
0.04458576
0.12546575
0.12571160
0.04921301
0.10151753
0.10491328
0.34646210
16
0.24253564
I8
0.22941573
-2.86674523
0.04912448
0.10688034
0.10763902
20
0.2 1821788
-89.92565918
0.05485482
0.09624152
0.10446023
I V . l . 4 . The p-algorithm 1.4. The
2 13
p-algorithm.
Another technique is based on the construction of interpolating continued fractions of the form
If we take for g;(z) the ith convergent
and take z = 00 as extrapolation point, then we can put for i even:
6; = lim 30’2
g i ( z ) = do
+ d:! + . . . + di
Now the sequence b i can be calculated using reciprocal differences. Hence we can set up the following scheme:
k =0,1,2,. .. Zk-1. k = 1,2,3,... tk,l = a k - ak-1 Xk - Z k - j __ tk,j = tk-l,j-2 + tk,j-l - tk-1,j-1 tk,O
= ak
Zk
k=j,j+l
where in fact, in the notation of chapter 111, t k $ 3' = p [ Z k - j ,
. . .,zk]
Then clearly for a even
.,
bi = ~ [ z o , . .z;]
The values
tk,j
= ti,i
are ordered as in the following table
,...
and j = 2 , 3
,...
I V . l . 4 . The p-algorithm
214
to ,o t l ,o
tl,l
t2 ,o
t2,l
t2,2
t 3 ,O
t3,l
t3,2
with t k , o = a k and t k , - i = 0 for k = 0 , 1 , 2 , . . .. if we use { z ; } ~ ~as N interpolation points for the algorithm of Bulirsch and Stoer and {Z;};~IN with
as interpolation points for the p-algorithm, then for even 4:
(see problem (3)). Remark t ha t the computation of t i , i takes a smaller effort than th a t of 8y). For more properties of the p-algorithm we refer to [4]. We illustrate the influence of the choice of the {z~},.~N with some numerical examples. In the tables 4.10. and 4.11. several { b i } i E m were constructed, respectively for {ai);,N
={ l / m } i E N
and {a;)ieN = {I/('
+l ) } a ~
All computations were performed in single precision arithmetic. In each table the first sequence {b;}i,m is constructed with xi = &,the second with z i = i, the third with Z; = i2 and the fourth with xi = e i . The choice of the 2; greatly influences the convergence behaviour of the resulting sequence {b;};,N. We see that we do not necessarily obtain better results when the 2 ; converge faster to infinity. More information on this matter can be found in [8].
215
I V . l . 4 . The p-algorithm
Table 4.10.
___ i
-t ,. b; = t i , ; 6. b; 2 ti,; 1 - *,I 2; = i 2 z’- -2; = i a; I -_ _ ~ -
6
b-=t., x; = e' 1
1,t
0
1.00000000
1.00000000
1.00000000
1.oooooooo
1.00000000
2
0.57735026
7.07813832
0.24 11 8 0 8 8
0.50412308
0.48505624
4
0.44721358
0.02822608
0.11828871
0.30880488
0.30424682
6
0.37786450
-0.00126043
0.07431557
0.22015348
0.34178661
8
0.33333334
-0.00124686
0.06307473
0.17057084
0.30704348
10
0.30151135
-0.004163Q2
0.05177462
0.1 4054570
0.28143203
12
0.27735010
0.04758168
0.12608825
0.26141584
14
0.25818889
0.04435266
0.1 1 4 5 1 1 1 3
0.24517243
16
0.24253564
0.00017822
0.04088086
0.10530381
0.23163427
18
0.2284 1 5 7 3
0.00026065
0.05084635
0.00724047
0.22012024
0.21821788
0.00000440
0.03507820
0.08777263
0.21018887
20
0.00010506 -0.00047068
Table 4.11. a; =
i __
b*.
ai
0
1.00000000
2
0.33333334
4
0.20000000
6
0.14285716
8
0 . 1 1111111
10
0.08080008
=z
1 i+l
t.. 1,'
xi = -
b; = t %. *. bi = t i , ; .2 2; = i -~ --xi = 9,
1.00000000
b . = t *. , i
xi = e'
1.00000000
1.00000000
0.00000006
0.25000003
0.2401564 1
0.00000162
0.00000003
0.08080016
0.1543827 8
0.00000038
0.00000002
0.04545450
0.31652598
-0 .OOOOOO 15
0.02703313
0.08417003
0.00000003
0.01783462
0.07815881
-2.41421294
-0.00000042 0.00000018
1.00000000
IV.1.5. The multivariate e-algorithm
216
In [49] a study was made of some nonlinear convergence accelerating methods. The +algorithm was found to be the most effective one in many cases. Therefore we now generalize the t-algorithm for the convergence acceleration of multidimensional tables. 1.5. The multivariate t-algorithm. can he considered as a table with single entry. For its convergence acceleration we constructed the univariate function
A sequence {ai},.,N
c 00
f(z) =
VUiX'
i=O
with ai = 0 for i
< 0 and calculated
diagonal Pad6 approximants evaluated a t
x = 1. Let us now first consider a table (aj,k
}j,ken
with double entry and with A = limj,k.+OOaj,k = limj-w(limk+OO ajk) = limk-+w(limj+m ajk). To accelerate its convergence we introduce
We will again construct a sequence ( b i } i E ~ that will converge faster t o A in some cases [ll].To this end we now calculate bivariate Pad6 approximants for f ( z , y) and evaluate them at ( 5 , y) = ( 1 , l ) . If we denote by Ut =
j+k==t
Qjk
= at,O
+ at-1,l + . . + al,L-1 + a0,L
and start the c-algorithm with th e partial Bums of f( 1 , l )
IV.1.5. The multivariate 6-algorithm
approximants ri,i(z,y) evaluated a t
then the bivariate diagonal l’adi: (0) (e, y) = (1,l) are given by t2i . Hence we put bi
217
(0) = e2i
Let US generalize the idea for a table with multiple entry define
{ ~ j ~ . , , j ~ } ~ ~ We , , . . , ~ , ~ ~
W
f ( z l , .. .,zk) =
jl
, .. . ,j,=O
Vajl.,.j, Zil.
..zc
with
It is easy to prove that
Again multivariate Pad6 approxinlants for f ( z l , . . ., zk) can be calculated and evaluated a t (z,, . . .,zk)= ( 1 , .. ., 1) via the t-algorithm. Sirice
where now
the
6f)
for f(z1,. . .,Zk) are given by
IV.1.5. The multivariate €-algorithm
218
We illustrate this type of convergence acceleration with a numerical example. Suppose one wants to calculate the integral of a function G ( z ~ , ..,. Zk) on a bounded closed domain D of Rk. Let D = [0,1] x . . . x [0, I] for the sake of simplicity. The table { a i l . . . i k } j ~ , , , , , j can k ~ ~ be obtained for instance by subdividing the interval [O,1] in the lthdirection (t! = 1 , . . ., k) into 2jl intervals of equal length h~ = 2-Je(jt = 0, 1 , 2 , . . .). Using the midpoint-rule one can then substitute approximations
to calculate the ail...j, . We restrict oursclves again t o the case of two variables. With hl = 2-j and hp = 2-k we get
The column c0(4( l = 0, 1, . . .) in the 6-table given by
was also used by Gens [21] to start the €-algorithm for the approximate calculation of multidimensional integrals by means of extrapolation methods. He preferred this method to six other methods because of its simplicity and general use of fewer integrand evaluations. For
we have
[1
1
lim
j,k-cce
ajk =
1
dzdy = 2 In 2 = 1.386294361119891...
In table 4.12. one can find the a i , slowly converging to the exact value of the integral because of the singularity of the integrand in the origin. The 6; = r,,;(l, 1) converge much faster. For the calculation of the b i we need the c~ ( l = 0 , . . ., 2t).
IV.1.5. The multivariate €-algorithm
Hence b; should be compared with an information, i.e. j -tk = 2i.
ajk
which uses the same amount of
Table 4.12.
d
aii
1 2 3 4
1.166666606667 1.269017019048 1.325743700744 1.355532404415
219
bi = t2; (0)
1.330294906166 1.396395820203 1.386872037696 1.386308917778
220
IV.2.1. Iterative methods brrsed on P o d i approzimation
$2. Nonlinear equations.
Suppose we want to find a root xu of the nonlinear equation
f (4 = 0
Here the function f may be real- or complex-valued. If f is now replaced by a local approximation then a zero of th at local approximation could be considered as an approximation for x*. Methods based on this reasoning are called direct. One could also consider the inverse function g of f in a neighbourhood of the origin, if it exists, and replace g by a local approximation. Then a n evaluation of this local approximation in 0 could be considered as an approximation for x* since g(0)
.=z *
Methods using this technique are called inverse. 2.1. Iterative methods baaed on Pad6 approzimation.
Let
2;
be an approximation for the root xu of ri(x)
and let
Pi
= -(x) cli
be the Pad6 approximant of order ( m ,n) for f in z i . Then the next approximation x;+1 is calculated such that Pi(Zi+l)
=0
In case p ; ( z ) is linear ( m = 1) the value zi+l is uniquely determined. It is clear that this is t o be preferred for the sake of simplicity. A well-known method obtained in this way is Newton’s method (rn = 1, n = 0) which can be derived as follows. The Taylor series expansion for f(z) at zi is given by f(5) = f ( 5 ; ) f f ’ ( Z i ) ( Z
- 5 i ) + f‘ I ( x i ) ( 5 - X i ) 2 2
Hence the Pad6 approximant of order (1,O) for f a t
xi
+. . *
(4.2.)
equals
and we obtain (4.3.)
IV.2.1. Iterative methods boaed on Pad4 approzimation
221
Another famous method is Halley’s method based on the use of the Pad4 a p proximant of order (1,l) for f at X i : (4.4.)
Since the iterative procedures (4.3.) and (4.4.) only use information in the point z; to calculate the next iteration point z;+1 they are called one-point. The order of convergence of these iterative methods is given in the next theorem.
Theorem 4.4. If { z ~ } ~converges ~ N t o a simple root z* of f and if r i ( z ) is normal for every i then
For the proof we refer to [52]. Similar results were given in [IS], [43] and 1441. So the order is a t least rn n 1 and it depends only on the sum of m and ra. Consequently the order is not changed when Pad6 approximants lying on a n ascending diagonal in the Pad4 table are used. Iterative methods resulting from the use of (m,n) Pad6 approximants with n > 0 can be interesting because the asymptotic error constant C* may be smaller than when n = 0 [43]. The iterative procedures (4.3) and (4.4) can also be derived as inverse methods (see problem (5)). Let us apply Newton’s and Halley’s method t o the solution of
+ +
The root Z* = 0. We use xo = 0.09 as initial point. The next iteration steps can be found in table 4.13. Computations were performed in double precision accuracy (56-digit binary arithmetic).
222
IV.8.1. Iterative methods baaed on P a d 4 approzimation
Table 4.13.
i
Newton
Halley
Zi
2;
0
O.QOOOOOO0D-01
1
0.80997588D-01
2
0.6559965QD-01
3
0.43022008D-01
4
~.18500311D-01
5
0.34212128D-02
6
0.11703452D-03
7
0.13697026D -06
8
0.18760851D-12
9
0.35197080D-24
10
0.00000000D+00
0.90000000D -0 1 -0.40647645D-02 0.35821646D-06 -0.24514986D-18 0.00000000D+00
It is obvious that a method based on the use of (m,n) Pad4 approximants for f with n > 0 gives better results here: the function f has a singularity at z = 0.1. Observe that in the Newton-iteration zg is a good initial point in the sense that from there on quadratic convergence is guaranteed: 1zi+1
- 2.1 = 1zi+l I N c ' J z ~- z * I 2 for i 2
o
with C’ = 10. For the Halley iteration we clearly have cubic convergence from X I on. The formulas (4.3.) and (4.4.)can also be generalized for the solution of a system of nonlinear equations
which we shall write as
IV.d.1. Iterative methods based on Pad6 approzimation
223
Newton's method can then be expressed as [45]
where F’(xr),. . . ,xr)) is the Jacobiau matrix of first partial derivatives evaluated at (zy),. .. ,zt))with
Let us now introduce the abbreviations Fi = F ( z f ) ,. . . ,xk(i))
F , ! = F ( z (4 , ,... ,zk (4) To generalize Halley's method we first rewrite (4.4.) as
Then for the solution of a system of equations it becomes [I41
224
IV.2.1. Iterative method8 based o n P a d 6 appton'mdion
where the division and the square are performed componentwise and F " ( x 1 , ... ,x k ) is the hypermatrix of second partial derivatives given by
~~
a 2f k a2f k ... a2f k axkaxs ... axlaxk
dxlax2
.
.
I
az:
which we have t o multiply twice with the vector --F,!-'Fi. This multiplication is performed as follows. The hypermatrix F " ( x 1 , . . . , Z k ) is a row of k matrices, each le x k. If we use the usual matrix-vector nlultiplication for each element in the row we obtain a2fi
c
a2fk
c
k
asl ax;Yi ... i= 1 axlax;
Yi
...
k
i= 1
jin) vk
In [14] is proved that the iterative procedure (4.6.) actually results from the use of multivariate Pad6 approximants of order ( 1 , l ) for the inverse operator of ~ ( x l. ., . , x k ) at ( x v ) , . . . , x g ) ).
IV.8.1. Iterative method8 baaed on Pod6 approzimation
225
To illustrate the use of the formulas (4.5.) and (4.6.) we shall now solve the nonlinear system
{ (-
f 1 (‘J
Y) = -=+I/ = e-=-y
f&,y)
- 0.1
=0
- 0.1 = 0
which has a simple root at 12.0.1)) = (2.302585092994O46.. 0.
.
As initial point we take ( d 0 ) , g ( ’ ) ) = (4.3,2.0). In table 4.14.one finds the consecutive iteration steps of Newton’s and Halley’s method. Again Halley’s method behaves much better than the polynomial method of Newton. Here the inverse operator G of the system of equations F has a singularity near to the origin and this singularity causes trouble if we get close to it. For
we can write -0.5(ln(0.1 0.5(ln(0.1
+ u) + ln(O.l + v)) + u) - ln(O.l + v))
With (do), y(O)) = (4.3,2.0) the value do)= f2(z(O),y(O)) is close t o -0.1 which is close t o the singularity of G. For the computation of th e Pad6 approximants involved in all these methods the €-algorithm can be used. Another iterative procedure for the solution of a system of nonlinear equations based on the e-algorithm but without the evaluation of derivatives can be found in [5]. Since it does not result from approximating the multivariate nonlinear problem by a multivariate rational function, we do not discuss it here.
a
J9 0
Ne ton
0.43000000D+01
H Iev
,A4 0.20000000D+01
N N
0.43000000D+01
0.20000000D+01 0.5705728dD-kOO
1
-0.22427305D+02
-0.24720886D+02
0.287Q8400D+01
2
-0.21927303D+02
-0.24228888D-k.02
0.22495475D+Ol
-0.52625816D-01
3
-0.21427303D+02
-0.23720888D-kO2
0.23018229D+Ol
-0.44Q01947D-02
4
-0.20827303D+02
-0.23228888D+02
0.23025841D+01
-0.57154737D-05
5
-0.20427303D+02
-0.2272Q888D-kO2
0.23025851D-l-01
-0.Q7680305D--11
6
-0.1 QQ27303DS-02
-0.22228888D-kO2
0.23025851D+Ol
-0.17438130D--16
7
-0.18427303D+02
-0.21720888D-bO2
8
-0.18Q27303D+02
-0.2122Q888D-tO2
Q
0.18427303D+02
-0.20728888D)+02
10
-0.17827303D)+02
-0.20229888D-kO2
11
-0.17427303D+02
-0.19720888D-kO2
12
-0.16927303D+02
-0.10228888D+02
13
-0.16427303D+02
-0.18728888D+02
14
-0.15Q27303D+02
-0.18228888D-I-02
15
-0.15427303D+02
-0.17729888D+02
16
-0.14927303D+02
-0.17229888D+02
17
-0.14427303D+02
-0.16729888D+02
18
-0.13027303D+02
-0.1622Q888D-tO2
18
-0.13427303D+O2
-O.l572Q888D+02
20
-0,12!327303D+02
-0.15228888D-tO2
rn
y(4
3 e
m
J
?:
;r a
R e
Iv.8.B. Iterative methods based on rational interpolation
227
3.2.Iterative methods baa e d o n rut i o n a1 int erpol at i o n. Let
with pi and 9, respectively of degree m and n, be such that in an approximation z i for the root z * of f
fr)(zi-j)
with m that
= f(')(Zi-j)
I!
= 0 , . . .,8j - 1
(4.7.)
-+ n + 1 = C'tz08 ~ Then . the next iteration step z;+] is computed such Pi(zi4-1)
=0
For the calculation of zi+l we now use information in more than one previous point. Hence such methods are called multipoint. Their order of convergence can be calculated as follows. Theorem 4.5.
If
(Z;};~Nconverges to a simple root z* of f and f(n+n+l)(z) continuous in a neighbourhood of z* with
with n
> 0 is
where f ( k ) ( z *= ) 0 if k < 0, then the order of the iterative method based on the use of r,(z) satisfying (4.7.) is the unique positive root of the polynomial .i+1
- 8ozj- B1zj-l
The proof can be found in 1521.
- .. .- 8j = 0
228
IV.2.2.Iterative methods baaed on rational interpolation
If we restrict ourselves to the case 8[
=8
l = 0 , . , ., j
then it is interesting t o note that t h e unique positive root of
!=O
increases with j but is bounded above by 8 + 1 [53 pp. 46-52]. As a conclusion we may say t h a t the use of large j is not recommendable. We give some examples. Take m = 1, n = 1, 8 = I and j = 2. Then z;+1 is given by
The order of this method is 1.84, which is already very close to 8 + 1 = 2. Take m = 1, n = 1, 80 = 2, 81 = 1 and j = 1 . Then x;+1 is given by
The order of this procedure is 2.41. The ease m = I , n = 0,8 = 1 and j = 1 reduces to the secant method with order I .82. Let us again calculate the root of
with initial points close t o the singularity in x = 0.1. The successive iteration steps computed in double precision (56-digit binary arithmetic) are shown in table 4.15.
IV.2.2. Iterative method8 based on rational interpolation
229
Table 4.15.
(4.8.)
(4.0.)
secant method
Xi
zi
Xi
d
0.80000000D-01
0
1 2 3
0.90OOOO0OD-01
0.80000000D-01
0.80000000D-01
0.85000000D-01
O.QO0OOOOOD-01
O.QOOO00OOD-01
-0.15847802D-03 0.15131452D-06 -0.51538368D-13
4
0.20630843D-24
5
0.00000000D+00
- 0.17324153D-03
0.71095917D-01
0 . 4 6 6 2 1186D-10
0.64701421D-01
-0.6139231 1 D-25
0.46636684D-01
0.00000000D+00
0.30206020D-01 0.14080364D-01
6
0.425123451)-02
7
0.59843882D-03
8
0.28 4 30 08 1 D -0 4
0
0.15223576D-06
10
0.38727362D--10
11
0.58956894D-16
12
0.22832450D-25
By means of the multivariate Newton-Pad4 approximants introduced in section 8.2. of chapter 111 the previous formulas can be generalized for the solution of systems of nonlinear equations. We use the same notations as in chapter 111 and as in the previous section. For each of the multivariate functions f,(q,.. .,zk) with j = 1 , . ..,k we choose
D = N = ((0,. . . , O ) , ( l , O , . . .,O ) , ( 0 , 1 , 0 , .. . , O ) , . .., (0,...,o, 1)) H={(2,0,
..., 0 ) , ( 0 , 2 , 0,...10),..., (0,..., 0,2)) 2 INk
c Nk
IV. 2.2. Iterative methods baaed on rational interpolation
230
Here the interpolationset N U H expresses interpolation conditions in the points
zp)). . ., ($',
..
.)
Remark that this set of interpolation points is constructed from only three successive iteration points. The numerator of
with possible coalescence of points, is then given by
NO ,...,O ( 21, . . ,zk) c02,00,. ..(00
coo, ...,00,02
where
Nl,O
,...,0 (21> .
1
* 3
c12,00,...,00
0
zk
. . . NO ,...,0,1(21, . . ... 0 ...
*
coo, ...,00,12
>
zk)
IV.B.8. Iterative methods based o n rational interpolation
231
The values e,,tl,...,,Ltk are multivariate divided differences with possible coalescence of points. Remark that this formula is only valid if the set H provides a systern of linearly independent equations. The next iterationstep (21 ( i + l ), . . ., zk (i+l)) is then constructed such that
p i k ( z y l ) ., .
.J
Z k( i + l ) ) = 0
For k = 1 and without coalescence of points this procedure coincides with the iterative method (4.8.). With k = 2 and without coalescence of points we obtain a bivariate generalixation of (4.8.). Let us use this technique t o solve the sgslem
{
e--l+Y e-Z-9
= 0.1 = 0.1
with initial points (3.2, -0.95), (3.4, -1.15) and (3.3, -1.00). The numerical results computed in double precision (56-digit binary arithmetic) are displayed in table 4.16. The simple root is (2.302585092994046.. . ,O.). In this way we can also derive a discretized Newton method in which the partial derivatives of the Jacobian matrix are approximated by difference quotients
N = H = ( ( 0 , ...,Q ) , ( l , O , ..., 0 ),,..J(Ol...lO,l)}
D
= ( ( 0 , . . .,n))
If we call this matrix of difference quotients AFi, then the next iterate is computed by meana of
As an example we take the same system of equations and the same but fewer initial pointa as above. The consecutive iteration steps computed in double preciaion (56-digit binary arithmetic) can now be found in table 4.17.
232
IV.2.2. Iterative method8 baaed on rational interpolation Table 4.16.
i
~
y")
~-
0.320000OOD+Ol
-0.Q50000OOD-l-00
0.34000000D-k01
-0.11500000Df01
0
0.33000000D-k01
-0.10000000D~01
1
0.25249070D-kOl
-0.22072875D-kOO
2
0.22618832D-tOl
3
0.231278OQD)+Ol
-0.101644QOD-01 -0.51269373D-03
0.41 9 7 l Q 4 4 D - 0 1
4
0.23030978D-kOl
5
0.23025801D-kOl
0.40675854D-05
6
0.28025851Df01
-0.2568682QD-08
7
0.23025851 D+Ol
-0.1277881 6D-13
8
0.23025851Df01
-0.1 1350Q32D-16
Table 4.17.
i 0 1
~~
0.34000000D-k01
-0.11500000D-k01
0.33000000D-b0 1
-O.lO0OOOOOD+Ol
-0.29618530D-kOO
0.21743833Df01 0.20884Q33D-kOl
2
0.32743183D-l-01
3
0 . 2 2 1 1 4 2 1 1 Di-0 1
-0.84011352D+01
4
0.3651533RDt-0 1
-0.72149651D-kOl
5
-0.17Q00083D-kO4
0.208541 11 D i - 0 4
divergence
IV.2.3. Iteratave method3 wing continued fractions
233
The rational method is again giving better results. Now the initial points are
such t h at u = f l ( z , y ) is close to -0.1 which is precisely a singularity of the inverse operator for the considered system of nonlinear equations. For a more stable variant of thc diacretized Newton method we refer to [24]. 2.3. Iterative method8 wing continued fractions.
If ri(z) is t h e rational interpolant of order (m, 1) for $(z), satisfying 1 ri(zw)) = --(&-L))
L = 0,.. . , m + 1
f
then ri(z) can be written in the form rj(z) = do
+ d l ( z - z(i-m-')
+ ldm(z
z(i-m-
-
) + . . . + dm-l(Z - z(i-m-1) 1)
). . .(. - z(i-2J)l
t
+I
2
-
1.. .(z- z ( i - 3 ) )
&-q
&+I
The coefficients d , ( j = 0 , 1 , . . ., m ) are divided differences while dm+l is an inverse difference. The root of ,!;(z) can be considered as an approximation for the root z* of f . SO z(i+l)
-2
(i-1)
- dm+l
This method can be compared with methods based on the use of rational interpolants of order (1, rn) for f(z).
2.4. The qd-algorithm. First we s t a t e an important analytical property of the qd-scheme. To this end we introduce the following nokition. For the function f(s) given by its Taylor series development f(.) = co c 1 2 c 2 22
+
we define the H ankel-det erminant s
+
+
234
IV.Z.4. The gd-algorithm
We now call the series f(z)ultimately k-normal if for every R with 0 5 n 5 k there exists an integer M,, such t h a t for m 2 M,,the determinant Ifm,,,is nonzero.
Theorem 4.8. Let the meromorphic function f(z) be given by its Taylor series development in the disk B(0,r) = {z E C 121 < r ) and let the poles wi of f in B(0, r ) be numbered such that
1
0
< Iw] L
Iwzl
I ... < r
each pole occuring as many times as indicated by its order. If the series f(z) is ultimately k-normal for some integer k > 0, then the qd-scheme associated with this series has the following properties: a) for each n such that 0 < n 5 k and such t h a t Iwn-~l < (tun( < Iw,,+l(, where wg = 0 and, if f has exactly k poles, wk+l = 00, we have Iim m-+m
b) for each n such t h a t 0
1
qLm) = Wn
and such that Iw,(
Iim
m-cc
< Iw,,+lI,
we have
eim) = O
The proof can be found in [29 pp. 612-6131. As a consequence of this theorem the qd-algorithm is an ingenious tool t o determine the zeros of polynomials and entire functions. For if p(z) is a polynomial then the zeros of p ( z ) are the poles of the meromorphic function f(z) = b(z). AIIYq-column corresponding to a simple pole of isolated modulus would tend to the reciprocal value of that pole. It would be flanked by e-columns that tend t o zero. If moreover f ( z ) is rational the last e-column would be zero. Unfortunately the qd-scheme as generated in section 3 of chapter TI is numerically unstable. Rounding errors play an important role due to the divisions by small e-values. However, the rhombus rules (2.6.) and (2.7.) defining the scheme may be rearranged as follows:
IV.2.4. The qd-algorithm
235
In this way quotients are formed only with the quantities q p ) which do not tend t o zero and the qd-scheme is generated row by row. This is called the progressive form of the qd-algorithm. The problem of getting started is solved when we also calculate q y t l ) and elk+’) for negative values of k, i.e. if we calculate the extended qd-scheme as described in section 3 of chapter TI. In [XY pp.641-6541 it is also shown how to deal with the case in which several poles have the same modulus. To conclude this section we illustrate the preceeding theorem with the values qim) and el") for f(z) =
sin x 0.1 - x
which is an ultimately 1-normal function. Clearly with wo = 0,
w1
= 0.1 and
1 lirn qIm) = 10 = Wl
m-w
lim
m-oo
elm)
=0
Table 4.18. 0.10000000000000E+02
e
0 . Q Q S333 33 3 3333 3 E 4-0 1
ey)
=
0 . 1 ~ ~ ~ 6 6 6 ~ 6 6 ~ 6 ~ 7 ~ - n i
o.iooooooooooonnrc-t-o2
el")
=
0.83472454091016E-05
0.10000008347245E~+02
t:
=
-n.s3172454ooioie~-o5
0.100000000000001.;+02
$1
=
-0.108743778QQi25E-08
O.QOOOOQQQQS0l26E+O1
ey)
=
0.10000000000000E+02
e?
O . ~ O ~ O O O O O O O O O O ~ E + ~t:?) ~ 0.10n00000000000E-t02
-0.16686866666687E-01
-
=
0.1Q8743776QQ125E-08 0.27600144302181E-12 -0.276001 44382181 E-12
IV.R.5. The generalized qd-algorithm
238
2.5. The generalized qd-algorithm. The function f(s)can also be given by its Newton series 00
f(z)=
C f[zo, -
7
~ i ] ~ i ( z )
i=O
where the ~ [ z o. ., . , q] are divided differences with possible coalescence of points and where
ll j=l (z - zj-1) i
Bi(s) =
If we define the generalized Hankel-determinants
with Rm,o= 1 then we call the Newton series ultimately k-normal if for every n with 0 5 n 5 k there exists an integer Adn such that for m 1 Mn the determinant Hm,nis non-zero.
Theorem 4.7. Let the sequence of interpolation points ( 5 0 , 5 1 , 5 2 , .. .>be asymptotic to the sequence { ZO, zl,.. . , z,, 20, z1, . . . ,z i , 2 0 , z1, . . . ,z j , . . .} in the sense that
lim
k-oo
sk(j+l)+i =
d = 0 , . . . ,j
z;
Let the function f(s) be meromorphic in
..
~ ( 2 0 , ., Z j , r )
= {z E C
I
- zo)(z - 21).
. .( 5 - zj)l 5 r)
and analytic in the sequence of points { s ~ } ~and ~ Nlet the poles w i of B(z0,. . . , z j , r ) be numbered such that for wi
we have
=
I(Wi
- 20)
0 < w1 Iwp
1 . .
(Wi
- 2j)l
s . .. < r
f in
IV.2.5. The generalized qd-algorithm
237
where the poles are counted with their multiplicities. If the Newton series is ultimately k-normal for some integer k > 0, then the generalized qd-scheme associated with this series has the following properties : a) for each n such that 0 < n 5 k and such that w,-~ < w, < w,+1 where wo = 0 and, if f has exactly k poles, wk+l = r , we have
b) for each ra such that 0 < n 5 k and such that w , < w,+1, we have
For the proof we refer to [9].
238
IV.3.1. The we of Pad4 approzirnants
53. Initial value problems. Consider the following first order ordinary differential equation: dv
- =j(z,y)
dz
for
2
E [a,bj
(4.10.)
with y(a) = yo. When we solve (4.10.) numerically, we do not look for an explicit formula giving y ( z ) as a function of z but we content ourselves with the knowledge of y(zi) a t several points zi in [a, 61. If we subdivide the interval [a,b ] , k
U Izi-1 ,z i ]
[a, 61 =
i= 1
where zi =a
+ ik
i = 0,. . ., k
with
b-a h=--
k
for k
>0
then we can calculate approximations yi+l for y ( z i + ~ )by constructing local approximations for the solution y ( z ) of (4.10.) a t the point 2 ; . We restrict ourselves now t o methods based on the use of nonlinear approximations.
3.1. The use of Pad6 approximantcr. Let us try the following technique. If 8 i ( z ) is the Pad6 approximant of a certain order for y(z) a t z i then we can put Yi+l
=ai(Zi+l)
which is an approximation for y ( z i + l ) . For the calculation of ai(z) we would need the Taylor series expansion of y(z) at z;, in other words
Since the exact value of y(z;) is not known itself, but only approximately by y i , this Taylor series development is not known and hence this technique cannot be applied. However, we can proceed as follows. Consider the power series
IV.3.1. The u8e of Pad6 approzirnants
239
Let r i ( z ) be the Pad6 approximant of order ( m , n )for this power series. If we put z = zi+l, in other words z - x i = h , we obtain 8' = 0,..., k - 1
yi+r = r i ( z i + I )
Hence we can write ui+l
= Y i + hg(zir Yip h)
i = 0,. . .,k - 1
(4.11.)
where g is determined by r i . Such a technique uses only the value of zi and y i t o determine yi+1. Consequently such methods are called one-step methods. Moreover (4.11.) is an explicit method for t he calculation of v i + l . It is called a method of order p if the Taylor series expansion for g(z,y,h) satisfies Y ( z i + l ) - u ( z i ) - hg(ziJ Y ( z i ) ,h) = 0(hp+')
Clearly (4.11.) is a method of order (m + n) if r i ( z ) is a normal Pad6 approximant. The convergence of (4.11.) follows if g(z, y , h) satisfies the conditions of the following classical theorem [30 p. 711.
Theorem 4.8. Let the function g(zyy , h) be continuous and let there exist a constant L such that
g('J YJ
= f(’,
Y)
is a necessary and sufficient condition for the convergence of the method (4.11.), meaning that for fixed z E [ a ,b ] .
2 40
IV.3.1. The we of Podi approzarnants
From the fact that t i is a Pad4 approximant it follows th a t the relation g(z,y, 0) = f(z,y) is always satisfied (see problem (8)). The case n = 0 results in the classical Taylor series method for the problem (4.10.), If we take rn = n = 1 we get (4.12.)
If y(z) is a rational function itself, then using (4.11.) we get the exact solution Yi+l = Y(z;+I)
at least theoretically, if the degrees of numerator and denominator of r i ( z ) are chosen in a n appropriate way. Techniques based on the use of Pad4 approximants can be interesting if we has a large negative real part consider stiff differential equations, i.e. if [ZO]. An example of such a problem is the equation dY = xy dz
(4.13.)
-
with Re(1) large and negative. Since the exact solution of (4.13.) is
we have lim y ( z ) = lim
z-00
2-00
2% =o
and we want our approximations y; to behave in the same way. Dahlquist [15] defined a method to be A-stable if it yields a numerical solution of (4.13.) with Re(1) < 0 which tends to zero as i -+ 00 for any fixed positive h. He also proved t hat there are no A-stable explicit linear one-step methods. Take for instance the method of Euler (rn = 1, n = 0): Y i + l = ~i
+ hf ( z i j y i )
Y i + l = (1
+ hX)Yi
We get
= (1
+ hX)'+1 yo
IV.3.2. The
t~t: of
rational anterpolanta
241
Clearly lim yi+l = yo Iim (I i-
i-+W
00
+ /A)'+' = o
only if
)I
+ hX( < 1
So for large negative X the steplength h has to be intolerably small before acceptable accuracy is obtained. lo practice h is so small t h a t round-off errors and computation time become critical. The problem is t o develop methods t h a t do not restrict the stepsize for slability reasons. If (4.11.) results from the use of t h e Pad6 approximant of order (rn,m), (m, rn 1) or (m,m + 2) then one gets an A-stable method [16]. This can be seen as follows. If f ( z , y) = Xy then r,(z) is the Pad6 approximant for the power series
+
Hence Y i + l = rm,rt(hA) ~i
with h = ~ i + - lz;, where rm+(x) is the Pad6 approximant of order (rn,n)for e 2 . A-stability now follows from the following theorem.
Theorem 4.9. If m = n or m = n - 1 or m = n - 2 then the Pad6 approximant of order (m,n) for ez satisfies
For the proof we refer to [3, 181. 3.2. The uae of rational interpolants.
It is clear t h a t if the interpolation conditions are spread over several points, then the computation of yi+l will need several xi-^ and yi-t(t = 0, 1, . . .). Such methods are called multistep methods. Let r i ( z ) be the rational Hermite iuterpolant of order (rn,n) satisfying
242
IV.8.2. The me of rational interpolants
where (j
+ I)(R + 1) = rn + n + 1
Here
Then an approximation for y(zi+l) can be computed by putting
This is a nonlinear explicit multistep method. A two-step formula is obtained for instance by putting j = 1, m = 2, n = 1, and 8 = 1 and using theorem 3.17.:
where f, = j(x;,y;) and fi-1 = f(z;-1, yi-1). We can also derive implicit methods which require an approximation the caIculation of y;+1 itself, by demanding
v;+l
for
where (j
+ 2)(8 + 1) - 1 = m + n +
For m = 1 = n, j = 0 and
8
1
= 1 we get the formula
(4.15 .)
IV.3.3. Predictor-corrector method8
243
For more information concerning such techniques we refer t o [36]. Remark th a t multistep methods are never selfstarting. Both explicit and implicit ( j + 1)- step methods are of the form Yi+l =
j
C
acyi-c
+ h g ( z i + l , . . ., z i - j l
yi+lJ..
-J
~ i - j lh)
L=O
and they have order p if the Taylor series expansion of g satisfies i
y(zi+t) -
C
aty(zi-t)
- hg(zi+l
J
.
‘1
zi-j] y ( z i + l ) j , ..,y(zi-j), h ) = 0(hp+l)
L= 0
Hence (4.14.) is a third-order method if the starting values are third-order and (4.15.) is second-order. When applied t o a stiff differential equation one should keep in mind that linear multistep methods are not A-stable if their order is greater than two. The following result is helpful if af&’yl is real and negative. We know th a t we can write for problem (4.13.) y(z) = v(2i-j)
e(~-~i-j)X
with Re(X) large and negative. €leiice it is interesting to take a closer look a t rational Hermitc interpolants for exp(z) in some real and negative interpolation points and also in 0. Theorem 4.10.
5 m and 8qi 5 n is such th a t ri( t ) (0) = 1 = exp(*)(O)for 0 5 t 5 rn + n - t2 with e 5 m
If ri(z) = pi(z)/qi(z) with api a)
b) ri(tt) = exp(&) for & then lri(z)l
< 1 if
m
< 0, 1 5 k’ 5 j
5 n and
with &
# &, 1 5 !# k 5 j
z is real and negative.
For the proof we refer to [32]. 3.3. Predict0 r-correct or methode. For a solution of the initial value problem (4.10.) we have y’(z) = f ( z , y(z)) for every z in [a, b]
244
IV.S.4. Numeracal results
If we integrate this equation on the interval [zi-i, z i + t ] with j , 8 2 0 we get PZitl
Now f can be replaced by an interpolating function, through the points j(zi-1 ,yj-l)), . . ., which is easily integrated. ( z 8 ,f(z;, yi)), If t! = 1 we get p r e d i c t o r - m e t h o d s because they are explicit. If 8 = 0 we get corrector-methods because the value of y; is needed for the computation of y;. These implicit formulas can be used t o update an estimate of yfz;)iteratively. W h e n f is replaced by an interpolating polynomial we get the well-known methods of Adarns-Bashforth ( j = 0 and !. = 1) and Adams-Moulton ( 3 = 1 and 6 = 0). When f is replaced by an intc3rpolating rational function we can get nonliriear formulas of predictor or corrector type. Sincc the rational interpolant must be integrated, it is not recommendable to choose rational functions with a denominator of high degree.
3.4. Nu me rical r es ult 8 . Let u s compare two Taylor series methods ( m= 2 and m = 3) with the explicit method (4.12.) for t h e solution of the equation y‘ = I
+ y2
for z
2o
Y(0) = 1
The theoretical solution is g = tg(z -t 7r/4). We take the steplength h = 0.05. As can be seen in table 4.19. the second order rational method gives even better results t h a n t h e third order Taylor series method, a fact which can be explained by the singularity of the solution y(z) at z = 7 . To illustrate A-stability we will compare Euler’s method ( m= 1,n = 0) with the formulas (4.12.))(4.14.) and (4.15.) for the equation y’ = -25y
for z E [0,I]
Y(0) = 1
The solution is known to be y = exp(-25z). For the results in table 4.20. we chose the steplength h = 0.1. As expected Eiiler’s solution blows up while formulas based on the use of a “diagonal” entry of the Pad6 table or the rational Hermite interpolation table for t h e exponential decay quite rapidly. Similar results would have been obtained if “superdiagonal” entries of the Pad6 table were used. Remark that formula
Table 4.19. exact solution
+7)
Taylor series
Pad6 approximant
Taylor series
m=3
i
zi
1
0.05
0.1 10536D+01
0.1 105OOD+Ol
0.110526D+01
0.110533D+01
2
0.10
0.122305D+01
O.l22210D+Ol
0.122284D-kOl
0.122299D)+Ol
3
0.15
0.135600D-kOl
0.135449D+Ol
0.135573D+01
0.135598D)+01
4
0.20
0.150850D+01
0.1 5 0 5 8 2 D t 0 1
0.150795D+01
0.150831D+Ol
5
0.25
0.1 6858oD+oi
0.168150D-kOl
O.l68500D+Ol
0.188547D+01
6
0.30
0.188577D$-01
O.l88886D+Ol
0.188462D-kOl
0.189522D-t-01
7
0.35
0.214875D-kOl
0.213884D-kOl
0.214811D+Ol
0.214883D-kOl
8
0.40
0.246496D-kOl
0.244751D-kOl
0.24626lD+Ol
0.246335D-kOl
9
0.45
0.286888D+01
0.28398OD+Ol
0.286543D-I-01
0.286584D-I-01
10
0.50
0.340822D+Ol
0.335737D-t-01
0.340298D-t-01
0.340248D-kO 1
11
0.55
0.4 16936D4-01
0.407388D-kOl
0.416097D+Ol
0.415703D+01
12
0.60
0.533186D-t-01
0.513307D+Ol
0.531720D+Ol
0.5301 31D-t-01
13
0.65
0.7 3 4 0 4 4 D t 0 1
0.685144D-t-01
0.731087D+01
0.724568D-kOl
14
0.70
0.1 16814D+02
0.1 00697D-kO2
0.1160 1QDi-02
0.112431D-kO2
15
0.75
0.282383D-kO2
0 .I 77676D-kO2
0.277486D-i-02
0.232131D-kO2
tg(zi
m=2
m=l=n
*CD e
5 e
Table 4.20.
exact solution exp(-25si)
i 1
0.1
0.820850D-01
2
0.2
0.6737851)-02
3
0.3
0.553084D-03
4
0.4
0.453QQQD-04
5
0.5
0.372665D-05
6
0.6
0.305902D-06
7
0.7
0.251 100D-07
8
0.8
0.2061 15D-08
Q
0.Q
0.16QlQOD-OB
10
1.0
0.13887QD--10
m=l=n explicit one-step
Euler -0.150OOOD+01 0.225000D+01 -0.337500D+01 0.506250D-l-01 -0.758375D+01 0.113Q06D+02 -0.17085QD)+02 0.25628OD-l-02 -0.384434D)+02 0.576650D4-02
1
-0.111111D+OO 0.123457D-01 -0.137174D-02 0.152416D-03 -0.189351 D-04 0.188168D-05 -0.200075D-06 0.2323061)-07 -0.258117D-08 0.286797D-00
rn = 2,n = I
explicit multisteD
rn=l=n
implicit one-step
0.820850D-01
0.123047D+00
0.840728D-01
0.151 407D-01
-0.542641D-01
0.186302D-02
-0.4Q474OD-kOO
0.22023QD-03
-0.252173D+00
0.282073D-04
0.314055D-kOO
0.347083D-05
0.102175D+O1
0.427077D-06
0.16Ql37D-kOO
0.525507D-07
-0.lQl5Q6D-I-00
0.646622D-08
-0.6Q8115D-kOO
0.7Q5651D-OQ
IV.3.5. Systems of first order ordinary differential equatiow
247
(4.14.) based on a "subdiagonal" rational approximation is surely not producing an A-stable method. To obtain the results of table 4.20. by means of (4.14.) a second starting value yl was necessary. We took yl = exp(-2.5) = y(z1). For (4.15.) the expression f(zi,y;) = -25y; was substituted and y;+l was solved
from the quadratic equation. All these schemes can be coupled to mesh refinement and the use of extrapolation methods. If an asymptotic error expansion of y(zi) in powers of h exists, then the convergence of the sequence of approximations for y(z;), witb z i fixed, obtained by letting the stepsize decrease, ran be accelerated by the use of techniques described in section 1 [23]. 3.5. Systems of first order ordinary differential equations.
Nonlinear techniques can also be used to solve a system of first order ordinary differential equations
dzj = f j ( X , 2 1 , 2 2 , . . . , Z t ) dX
where the values z,,! and the functions approaches are possible. If we introduce vectors
fj
j = 1, ..., k are given for j = I , . . .,k. Several
then one method is t o approximate the solution componentwise using similar techniques as in the preceding sections. So for instance (4.12.) becomes
Y ( z ; + ~N) Y;:+l where
and
= Y;
Y;) ] + h [*F(zi,2F2(Zi, Y;)- hF'(xil yi)
(4.16.)
248
IV.3.5. Systems of first order ordinary differential equation8
and the addition and multiplication of vectors is performed componentwise. In other words (4.16.) is equivalent with
For more information on such techniques we refer to [56,35, 381. Another approach is not based on componentwise approximation of the solution vector Y ( s )but is more vectorial in nature. Examples of such methods and a discussion of their properties is given in [57, 7 , 27, 131. The nonlinear techniques introduced here c a n also be used to solve higher order ordinary differential equations and boundary value problems because these can be rewritten as systems of first order ordinary differential equations. Again same of the nonlinear techniques prove to be especially useful if we are dealing with stiff problems. A system of differential equations
dY -= F ( z , Y ) dx Y ( a ) = Yo with
is called stiff if the matrix
IV.3.5. System8 of first order ordinary differential equations
249
has eigenvalues with small and large negative real part.
{
Consider for example
+
= 9 9 8 ( ~x )~ 199822 (z)
Zl(0j = 1
- - - -99921(2) - 19992z(x)
22(0)
=0
The solution is
so that again
lim z I ( x ) = 0 = lim
z 4 m
DO’%
z2(x)
where both z1 and a2 contain fast and slow decaying components. For a discussion of stiff problems we refer to [20 pp.209-2221.
IV.4. Numerical integration
250
§4. Numerical integration.
s,
b
Consider I = j ( z ) d z . Many methods to calculate approximate values for I are based on replacing j by a function which can easily be integrated. The classical Newton-Cotes formulas are obtained in this way: j is replaced by an interpolating polynomial and hence I is approximated by a linear combination of function values. In some cases the values of the derivatives of j ( x ) are also taken into consideration and then linear combinations of the values of f ( z ) and its derivatives a t certain points are formed t o approximate the value I of the integral. This is for instance the case if polynomial Hermite interpolation is used. In many cases the linear methods for approximating I give good results. There are however situations, for example if f has singularities, for which linear methods are unsatisfactory. So one could t r y t o replace f by a rational function r and consider
[
r(x)dz
as an approximation for I . But rational functions are not that easy to integrate unless the poles of r are known and the partial fraction decomposition can be formed. Hence we use another technique. Let us put
Then
1 = y(b) If f is Riemann integrable on [a,b ] , then y is continuous on [a, b ] . If f is continuous on [a, 61, then y is differentiable on [a,b] with y’(z) = j ( z ) and y(a) = 0
So I can be considered as the solution of a n initial value problem and hence the techniques from the previous section can be used. We group them in different categories.
I V . 4 . l . Methoda w i n g Pad6 approximantd
251
4.1. Methods using Pad6 approzirnants. Let us partition the interval [a, 61 with steplength h = ( b - a ) / k and write Z; = a
+ ih
a = 0, ...,k
and (4.17.)
s,”’
where y ; approximates y ( z i ) = f(t)dt. If r; is the Pad6 approximant of order ( m ,n ) for t l , ; ( h )then we can put
and consider y k as an approximation for I . In this way Yi+l
which means
= yi
+ hg(zi, h)
i = 0, ...,k - 1
(4.18.)
l;’’
f ( t ) d t N hg(z;,h)
If m = 1 = ra we can easily read from (4.12.) t h a t (4.18.) results in (4.19.) Formulas like (4.18.) use derivatives of f(z) and are nonlinear if n > 0. From the previous section we know that (4.18.) is exact, in other words that yk = I , if y ( z ) is a rational function with numerator of degree m and denominator of degree ra. For n = 0 formula (4.18.) is exact if y ( z ) is a polynomial of degree m, i.e. f ( z ) is a polynomial of degree m - 1. The obtained integration rule is then said t o be of order m - 1. The convergence of formula (4.18.) is described in the following theorem which is only a reformulation of theorem 4.8.
252
IV.4.2. Method8 wing rational interpolanta
Theorem 4.11. Let y; be defined by (4.18.). Then lim h-0
yi(h) = y(z) for fixed z E [a, b]
====i(h)
if and only if g(z,O) = f(z). Instead of (4.17.) one can also write t l , i ( h )= T/i
+ lata,i(h)
with
and compute Pad6 approximants r , of a certain order for tz,;(h).Let us take
m = 1 = n . If we define
Yi+i
= yi
+ hr;(h)
then we get
4.2. Method8 using rational interpolanta.
If we again proceed as in the section on initial value problems we can construct nonlinear methods using information in more than one point. Since these methods are not self starting but need more than one starting value their use is somewhat limited and rather unpractical. An example of such a procedure is the following. Let ti(.) be the rational Hermite interpolant of order ( 2 , l ) satisfying
and let
IV.4.3. Methods without the evaluation of derivatives
253
Then we know from formula (4.14.) t h a t
with f; = f(zi) and f i - 1 = f ( z i - I ) . Often rational interpolants are preferred t o Pad6 approximants for the solution of numerical problems because the use of derivatives of f is avoided. As mentioned, a drawback here is the necessity of more starting values. Another way t o eliminate the use of derivatives, now without the need of more starting values, is the following.
4.3. Methods without the evaluation of derivatives. One can replace derivatives of j ( z ) in formula (4.18.) by linear combinations of function values of f ( z ) without disturbing the order of the integration rule. To illustrate this procedure we consider the case m = 1 = n. Then 1591
We will compute constants a,p and 7 such t h a t for
we have gfZ,
h) - t ( s ,h) = O(hrn+%) = O(h2)
For ~ ; + i= ~i
+ h t ( z ; ,h)
this would imply y(si+l) - Yi+l = ~(hm+m+l)= O ( P )
Condition (4.20.) is satisfied when a+B=2 B7 = -1
(4.20.)
254
IV.4.4. Numerical resulta for singular integranda
In other words, for 7 # 0, 27 a=---
p = -1
+1 7
-.
7
so
For 7 = 1 we get the integration rule
In this way we approximate (4.21.)
4.4. Numerical results for singular integranda We will now especially be interested in integrands regular in [a,b] but with a singularity close t o the interval of integration and on the other hand in integrands singular in Q or b . The problem of integrating a function with several singularities within [a, 61 can always be reduced t o the computation of a sum of integrals with endpoint singularities. If f ( z ) is singular in 6 then the value of the integral is defined by
and is assumed t o exist. We shall compare formula (4.19.) with Simpson’s rule ( m = 2, n = 0 ) and with a (2k)-point Gaussian quadrature rule that isolates the singularity in the weight function. If f(t) can be written as w ( t ) h ( t ) where w ( t ) contains the endpoint singularity of f(t) and h(t) is regular then the approximation
I
wlh(tl) + . . .
+ 402kh(t2k)
does not involve function evaluations in singular points. We use a (2t)-point formula because on [zi,2;+1] for i = 0 , . . ., k - 1 both (4.19.) and Simpson’s rule
IV.4.4. Numerical reaulta for aingular integranda
255
need two function evaluations. Since / is singular in b = zk, we take f ( z k ) = 0 in Simpson's rule which means that the singularity is ignored. In (4.19.) the singularity of f is no problem since f and f’ are only evaluated in 2 0 , . . ., xk--1 [58]. Our first numerical example is
d t = 3.04964677.. with a singularity in
t = In 3 = I .09861228. . . and the second example is
+
2(1_ _ ~ t)sint _ _ _ _ _cost d t = 2
We shall also compare the different integration rules for the calculation of e'dt = 1.71828182. . .
which has a smooth integrand. Because of the second example the weight function w ( t ) in the Gaussian quadrature rule was taken t o be W ( t )=
1
J1-t
.___
All the computations were performed in double precision arithmetic and the double precision values for the weights w i (i = 1 , . . .,2k) and the nodes ti (d = 1 , . . .,2k) were taken from [I pp. 916-9191. Remark t ha t the nonlinear techniques behave better than the linear techniques in case of singular integrands. However, for smooth integrands such as in table 4.23., the classical linear methods give better results th a n the nonlinear techniques. Also, if the singularity can be isolated in the weight function such as in table 4.22., Gaussian quadrature rules are very accurate. In general, little accuracy is gained by using nonlinear techniques if other methods are available for the type of integrand considered [51].
256
IV.4.4. Numerical reeults f o r singular integranda
As in the previous section all these schemes can be coupled to mesh refinement and extrapolation.
Table 4.21. 1
I=J
et
o (3 - et)2
k = 16
dt = 3.04964677.. .
Gaussian
Simpson's
formula
quadrature
rule
(4.19 .)
3.1734660 3.0773081 3.0564777
3.28067 65 3.0841993 3.0531052
3.1006213 3.0573798 3.0510553
Table 4.22.
I = J 2(1 - t ) sin t -tc o s t dt = 2 0 G - j
k=4 k=8 k = 16
Gaussian
Simpson's
formula
quadrature
rule
(4.19.)
2.0000000 2.0000000 2.0000000
1.7504074 1 A267581 1.8786406
1.9033557 1.8797832 1.go48831
Table 4.23. 1
I = 1etdt = 1.71828182.. . 0
k=4 k=8 k = 18
Gaussian
Simpson's
formula
quadrature
rule
(4.18.)
1.7205746 1.7204038 1.7188196
1.7182842 1.7182820 1.7182818
1.7284991 1.7206677 1.7188592
IV.5. Partial diflerential equatiom
257
§5. Partial differential equations. Nonlinear techniques are not frequently used for the solution of partial differential equations. We will describe here a method based on the use of Pad6 approximants t o solve the heat conduction equation which is a linear problem. For other illustrations we refer to [54, 48, 17, 22, 281. All these techniques first discretize the problem so that thr: original partial differential equation is replaced by a system of equations which is nonlinear if the partial differential equation is. Another type of techniques which we do not consider here are methods which do not discretize the original problem but solve it iteratively by means of a procedure in which subsequent iteration steps are differentiable functions [13]. Linear techniques of this type arc recommended for linear problems and nonlinear techniques can be used for nonlinear partial differential equations. Let us now concentrate on the heat conduction equation. Suppose we want t o find a solution u ( z , t )of the linear problem
a+, t ) - a%+, t ) ~
at
.~
82
~
a
t>O
(4.22 .)
with u(2,O) = u ( 2 )
u(0,t ) = a
4 4 t)= P The domain [a, 61 with
x
[0,oo) is replaced by a rectangular grid of points (z;, ti)
zi = a + ih ti = j A t
i = 0,.. .,k + 1
j = 0, 1 , 2 , . . .
where
and At is the discretization step for the time variable. We first deal with the discretization in the space variable z and introduce u;(t) = u ( z ; , t ) for
t20
Then using central differences, (4.22.) can for instance be reduced t o
258
I V.5. Part id digere nt ial
equ at ion8
for t > 0 and i = 1, . . . , k. In general we c a n write
(4.23.)
where A is a real symmetric positive definite k x k matrix and depends on the chosen approximation for the operator d 2 / d s 2 . If we introduce the notations
then the exact solution of (4.23.) is U ( t ) = eWtAV where e-tA is defined by
Using the discretization in the time variable t we can also write (4.24.)
I V.5. Partial d ifler e nt ial equ at i o ns
259
If rm,.,(t)= p ( t ) / q ( t ) is the Pad6 approximant of order (n,n) for e-t then (4.24.) can be approximated by
U ( t + A t ) = [q(At.A)]-' [P(At.A)] U ( t )
(4.25 .)
Varga proved t ha t for n 2 m this is an unconditionally stable method [ 5 5 ] , meaning t ha t initial rounding errors remain within reasonable bounds as the computation proceeds independent of the stepsize used. If rn = 1 and n = 0 then (4.25.) means
U(t + At) = ( I
- AtA)U(t)
or equivalently if A = (&.!)kXk
ui(ti+l) = u i ( t j ) - ~t
k
C aitut(tj)
.!=I
which is the well-known explicit method to solve (4.22.). The solution a t level t = ti+l is determined from the solution at t = t , , For nt = 1 = n we obtain from (4.25.)
or equivalently
which is the method of Crank-Nicholson [50]. The operator a2/ax2is replaced by the mean of a n approximation for the partial second derivative a t level t = and the same approximation at t = t i .
IV.6.1. Kernel8 of finite rank
260 56. Integral equations.
As in the previous section we shall discuss linear equations for which the use of nonlinear methods is recommendable. Those interested in nonlinear integral equations are referred to [lo, 121 where methods are indicated for their solution. If the integral equation is rewritten as a differential equation then techniques developed for the solution of initial value problems can also be used. We restrict ourselves t o the discussion of an inhomogeneous Fredholm integral equation of the second kind (the unknown function f appears once outside the integral sign and once behind it):
f ( 4-
J'." K ( z ,y)f(v)dy
= g(z)
xE
1% bl
(4.28.)
Here the kernel K ( z , y ) and the inhomogeneous right hand side g(z) are given real-valued continuous functions. Fredholm equations reduce to Volterra integral equations if the kernel K ( z , y ) vanishes for g > z which produces a variable integration limit.
6.1.K e r n e l 8 of finite rank. Formally the solution of (4.26.) can be written as a series. P u t
and
If we define
and
then (4.26.) reduces to
IV.6.1. Kernel8 of finite rank
261
The series (4.27.) which is a power series in X, is called the Neumann series of the equation (4.26.). Convergence of the Neurnann series for certain values of X depends on the properties of the kernel K ( z ,y). If K ( z ,y ) is bounded by IWz,y)l
for (z,Y) E
b,b]
X [a, b ]
then clearly the series (4.27.) converges uniformly t o f(s)in [a, b] if 1
1x1 < ---___ M ( b - a) If Pad6 approximants in the variable X are constructed for (4.27.) then they may have a larger convergence region than the series itself. Especially interesting is the case t ha t the kernel is degenerate, in other words k
~ ( zY), =
C Xi(Z)Yi(y)
i= 1
with ( X i } and {Y;} each linearly independent sets of functions. Such a kernel is also said to be of finite rank k. Let us try to determine f(z) in this case. If we put
then (4.26.) can be written as (4.28.) Multiply (4.28.) by Yj and integrate to get
or equivalently
IV.6.2. Completely continuoua kernela
262
(4.29 .)
which is the determinant of the coefficient matrix of system (4.29.), then D(X) is a polynomial in X of degree at most k. In case D(X) # 0 the solution of (4.29.) is given by
with
Dij
the minor of the ( i , j ) t h element in D(X). SO (4.28.) can be written as
which is a rational function in X of degree a t most k both in numerator and denominator. Since the series development of (4.30.) coincides with the Neumann series (4.27.) we know that the solution f(z) is equal to the Pad6 approximant of order (k,k) for the Neumann series. So in this case the Pad6 approximant is the exact sum of the series because the sum is a rational function [2 p. 1761. 6.8. Completely continuoua kernels.
The equation (4.26.) can be rewritten
EUI
(I-XK)f=g where the linear operators I and K are defined by
(4.31.)
263
IV.6.2. Completely continuoua kernela
We shall consider square-integrable functions f with
Suppose now t ha t { f;}iEN is a bounded converging sequence of functions. Then the operator K is said t o be completely continuous if for all bounded { f ; ) , ~ the sequence { K fi}iENcontains s subsequence converging t o some function h ( s ) with
- h(2)l2ds -+ 0 as i -+
llKfi - hll = \i[lKfi(Z)
00
A basic property of completely continuous transformations is th a t they can be uniformly approximated by transformations of finite rank. Thus there is a n infinite sequence { K ; } ; Lof~ kernels of finite rank i such th a t
where lim
i-+w
t;
=0
When a completely continuous kernel K is replaced by Ki then the solution f; of
( I - Mi)!; = g is given by the Pad4 approximant ri,i of order (i,;), as explained in the previous section. In case K is completely contiuuous the exact solution f of (4.31.)is a meromorphic function of X [31 p. 311 arid
lim f, = f
i-
00
2 64
Iv.6.2. Completely continuous kernels
in any compact set of the X-plane except at the finite number of X-poles of f [8].Since
we consequently have
lim
ri,i
1-m
=f
in any compact set of the X-plane except at the finite number of X-poles of f.Thus the solution f is the limit of a sequence of diagonal Pad6 approximants. This result is not very useful since each approximant in the sequence is derived from a different Neumann series, with kernel K ; and not with kernel K . However a similar result exists for the sequence r i , i derived from the Neumann series (4.32.) i= 1
where the operator K' is defined by
Theorem 4.12. If ri,i is the Pad4 approximant of order (d, i) t o the Neumann series (4.32.) of the integral equation (4.28.) with completely continuous kernel K then the solution f is given by lim
i-oc
ri,i
=f
in any compact set of the X-plane except a t the finite number of poles of f and at limit points of poles of ri,i as i -+ 00. More information on this subject can be found in [2 pp. 178-1821.
I V . Problems
265
Problems. Show t ha t with
EF’= an for n 2 0,
which coincides with Aitken’s A2-process to accelerate the convergence of a sequence. Give a n algorithm similar to the one in section 1.4. for the calculation of t,,i, but now based on the use of inverse differences instead of reciprocal differences, Show that if the algorithm of Bulirsch and Stoer is used with the interpolation points zi (lim,+m s; = 0) and the p-algorithm is used with the interpolation points z: = 1,s; (1imi4- 2;. = oo),then t i ; = B~( 0 ) . Compare the amount of additions and multiplications performed in the c-algorithm and the qd-algorithm when used for convergence acceleration. Derive the formulas (4.3.) and (4.4.) using inverse interpolation instead of direct interpolation. Derive the formulas (4.8.) and (4.9.) based on the use of rational interpolants. in section 2.3. for successive values of Organize the computation of d,+l d such t ha t a mininum number of operations is involved. Prove that one-step explicit methods for the solution of initial value problems based on the use of Pad6 approximants are convergent if g(z, g, h) given by (4.11.) is continuous and satisfies
IV. Problem8
266
(9)
Check the formulas (4.12.) and (4.14.).
(10)
Write down formula (4.16.) for the solution of
- 99821 (2)+ 199822( Z) - - - -9992+) (11)
- 199922(.)
ZI(0) = 1 ZZ(0) = 0
Construct a nonlinear numerical integration rule based on the use of Pad4 approximants of order ( 2 , l ) . Afterwards eliminate the use of derivatives as explained in section 4.3.
I V . Remnrka
267
Remarks. Nonlinear methods can also be used for the solution of other numerical problems. We refer for instance to [40] where the solution of linear systems of equations is treated, t o [19] for analytic continuation, to [34] for numerical differentiation and to (421 for the inversion of Laplace transforms.
An important link with the theory of numerical linear algebra is through QR-factorization. Rutishauser [47] proved th a t the determination of the eigenvalues of a square malrix A can be reduced to the determination of the poles of a rational functiom .f built from th e given matrix. In this way decomposition techniques for A to compute its eigenvalues are related t o the qd-algorithm when used to compute poles of meromorphic functions. Univariate and multivariate continued fractions and rational functions are also often used t o approximate given functions. For univariate examples we refer t o (33, 421. The bivariate Beta function is a popular multivariate example because it has numerous singularities in a quite regular pattern. For numerical results we refer to [89, 26, 121.
As a conclusion we can say t h a t every linear method has its nonlinear analogue. In case linear methods are inaccurate or divergent, it is recommendable t o use a similar nonlinear technique. The price we have to pay for the ability of the nonlinear method to cope with the singularities is the programming difficulty to avoid division by small numbers.
2 68 References.
1
1) Abramowite M. and Steguo Handbook o Mathematical functions. Dover publications, New York, 1968.
[
21 Baker G. and Gamrnel J. The Pad6 approximant in theoretical physics. Academic Press, New York, 1970.
[
31 33aker G. and Graves-Morris P. Pad4 approximants: basic theory. Encyclopedia of Mathematics and its applications: vol 13. AddisonWesley, Reading, 1981.
\
41 Brezinski C. Algorithnies d’accklkration de la convergence. Editions Tecbnip, Paris, 1978.
[
51 Rrezinski C. Application de 1’6-algorithme B la r6solution des syst6mes non linhaires. C. R. Acad. Sci. Paris Sbr. A 271, 1970, 1174-1177.
[
61 Bulirsch R . and Stoer J. Fehlerabschatzungen und Extrapolation mit rationalen Funktionen bei Verfahren vom Richardson-Typus. Numer. Math. 6, 1964, 413-427.
[
On the stability of rational Runge 71 Calvo M. and Mar Quemada M. Kutta methods. J. Comput. Appl. Math. 8, 1982, 289-292.
[
81 Chisholm J. Solution of linear integral equations using Pad6 Approximants. J. Math. Phys. 4, 1963, 1506-1510.
[
91 Claessens G. Convergence of multipoint Pad6 approximants. Report 77-26, University of Antwerp, 1977.
[ lo] Clarysse 7 .
Rational predictor-corrector methods for nonlinear Volterra integral equations of the second kind. In [60], 278-294.
[ 111 Cuyt A.
Accelerating the convergence of a table with multiple entry. Numer. Math. 41, 1983, 281-286.
[ 121 Cuyt A.
Pad6 approximants for operators: theory and applications. Lecture Notes in Mathematics 1065, Springer Verlag, Berlin, 1984.
[ 131 Cuyt A.
Pad6 approximants in operator theory for the solution of nonlinear differential and integral equations. Comput. Math. Appl. 6, 1982, 445-466.
Iv. Reference8
269
Abstract Pad6 approximants for 141 Cuyt A. and Van der Cruyssen 1’. the solution of a system of nonlinear equations. Comput. Math. Appl. 9, 1983, 617-624.
[ 151 Dahiquist G.
A special stability problem for linear multistep methods. BIT 3, 1963, 27-43.
[ 161 E61e B.
A-stable methods and Pad4 approximations t o the exponential. SIAM J. Math. Anal. 4, 1973, 671-680.
[ 171 Fairweather G.
A note on the eficient implementation of certain Pad6 methods for linear parabolic problems. BIT 18, 1978, 100-108.
[ 181 Frame J .
The solution of equations by continued fractions. Amer. Math. Monthly 60, 1953, 293-305.
[ 191 Gammei J.
Continuation of functions beyond natural boundaries. Rocky Mountain J. Math. 4, 1974, 203-206.
[ 201 Gear C .
Numerical initial value problems in ordinary differential equations. Prentice-Hall Inc., New Yersey, 1971.
[ 211 Genz A .
The approximate calculation of multidimensional integrals using extrapolation methods. Ph. 1). in Appl. Math., University of Kent, 1975.
[ 221 Gerber P. and Miranker W .
Nonlinear difference schemes for linear partial differential equations. Computing 11, 1973, 197-212.
[ 231 Gragg W.
On extrapolation algorithms for ordinary initial value problems. SLAM J. Numer. Anal. 2, 1965, 384-403.
[ 241 Gragg W. and Stewart G .
A stable variant of the secant method for solving nonlinear equations. SIAM J. Numer. Anal. 13, 1976, 889-903.
[ 251 Graves-Morris P .
Pad6 approximants and their applications. Academic
Press, London, 1973.
[ 261 Graves-Morris P. , Hughes Jones R . and Makiuson G . The calculation of some rational approximants in two variables. J. Inst. Math. Appl. 13, 1974, 311-320.
[ 271 Hairer E.
Nonlinear stability of RAT, an explicit rational Runge-Kutta method. BIT 19, 1979, 540-542.
IV. References
270
1
Pad6 approximants, fractional step methods 281 Hall C. and Porsching T. and Navier-Stokes discretizations. SLAM J . Numer. Anal. 17, 1980, 840851.
[ 291 Henrici P.
Applied and computational complex analysis: vol. 1. John Wiley, New York, 1974.
[ 301 Henrici P.
Discrete variable methods in ordinary differential equations. .John Wiley, New York, 1962.
[ 311 Hoheisef G.
Integral equations. Ungar, New York, 1968.
[ 32) Iserles A.
On the generalized Pad6 approximations to the exponential function. SLAM J. Math. Anal. 16, 1979, 631-636.
[ 33J Kogbethnts E .
Generation of elementary functions. In [46], 7-35.
[ 341 Kopal Z.
Operational methods in numerical analysis based on rational approximation. In 1371, 25-43.
[ 351 Lam bert J.
Computational methods in ordinary differential equations. John Wiley, London, 1973.
A method for the numerical solution of [ 361 Larnbert J. and S6aw B. g' = f ( z , y) based on a self-adjusting non-polynomial interpolant. Math. Comp. 20, 1966, 11-20.
[ 371 Langer R.
On numerical approximation. University of Wisconsin Press, Madison, 1959.
[ 381 Lee D.and Preiser S.
A class of nonlinear multistep A-stable numerical methods for solving stiff differential equations. Internat. J. Comput. Math. 4, 1978, 43-52.
I 391 Levin D.
On accelerating the convergence of infinite double series and integrals. Math. Comp. 35, 1980, 1331-1345.
[ 401 Lindskog G.
The continued fraction methods for the solution of systems of linear equations. BIT 22, 1982, 519-527.
[ 411 Longman I.
Use of Pad4 table for approximate Laplace transform inversion. In [25], 131-134.
The special functions and their approximations. Academic 421 Luke Y . Press, New York, 1969.
IV. References
271
[ 431 Merz G .
PadCsche Naherungsbruche und Iterationsverfahren hoherer Ordnung. Computing 3, 1968, 165-183.
[ 441 Nourein M.
Root determination by use of Pad6 approximants. BIT 16, 1976, 291- 297.
[ 451 Ortega J. and RheinboIdt W.
Iterative solution of nonlinear equations in several variables. Academic Press, New York, 1970.
[ 461 Ralston A. and Wilf S.
Mathematical methods for digital computers. John Wiley, New York, 1960.
Der Quotienten-Differenzen AIgorithmus. [ 471 Rutishauser H. Math. Phys. 5, 1954, 233-251.
z. h g e w .
On time-discretiEations for linear time-dependent partial differential equations. University of Manchester, Numer. Anal. Report 5, 1974.
[ 481 Siemieniuch J . and Gladwell I .
Numerical comparison of nonlinear convergence accelerators. Math. Comp. 38, 1982, 481-499.
[ 491 Smith D. and Ford W. [ 501 Smith G.
Numerical solution of partial differential equations. Oxford University Press, London, 1975.
[ 511 Squire W.
In defense of linear quadrature rules. Comput. Math. Appl. 7, 1981, 147-149.
[ 521 Tornheim L.
Convergence of multipoint iterative metho&. J . Assoc. Comput. Mach. 11, 1964, 210-220.
[ 531 n a u b J .
Iterative methods for the solution of equations. Prentice-Hall Inc., New York, 1964.
[ 541 Varga R.
Matrix iterative analysis. Prentice-Hall Inc., Englewood-
Cliffs, 1962.
[ 551 Varga R .
On high order stable implicit methods for solving parabolic partial differential equations. J . Math. Phys. 40, 1961, 220-231.
[ 561 Wambecq A.
Nonlinear methods in solving ordinary differential equations. J. Comput. Appl. Math. l , 1976, 27-33.
Iv. References
272
[ 571 Wambecq A.
Rational Runge-Kutta methods for solving systems of ordinary differential equations. Computing 20, 1978, 333-342.
[ 581 Werner H . and Wuytack L.
Nonlinear quadrature rules in the presence of a singularity. Comput. Math. Appl. 4, 1978, 237-245.
[ 591 Wuytack L.
Numerical integration by using nonlinear techniques. J. Comput. Appl. Math. 1, 1975, 267-272.
[ 601 Wuytack L.
Pad6 approximation and its applications. Lecture Notes in Mathematics 765, Springer, Berlin, 1979.
1
611 Wynn P . Acceleration techniques for iterated vector and matrix problems. Math. Comp. 16, 1962, 301-322.
[ 621 Wynn P. 175-195.
Singular rules for certain non-linear algorithms.
BIT 3, 1963,
.
CBAPTER N:Applications
5 1. Convergence 1.1. 1.2. 1.3. 1.4. 1.5.
The The The The The
acceleration . . . . . . . . . . . . . . . . . . . . univariate t-algorithm . . . . . . . . . . . . . . . qd-algorithm . . . . . . . . . . . . . . . . . . . . algorithm of Bulirsch-Stoer . . . . . . . . . . . . . p-algorithm . . . . . . . . . . . . . . . . . . . . . multivariate €-algorithm . . . . . . . . . . . . . .
197
. 197 205
. 209 213
. 216
$2. Nonlinear equations . . . . . . . . . . . . . . . . . 2.1. Iterative methods based on Pad4 approximation . 2.2. Iterative methods based on rational interpolation 2.3. Iterative methods using continued fractions . . . 2.4. The qd-algorithm . . . . . . . . . . . . . . . 2.5. The generalized qd-algorithm . . . . . . . . .
. . . . . 220 . . . . . . 220 . . . .
. . . . . 227 . . . . . 233 . . . . 233 . . . . . 236
53. Initial value problems . . . . . . . . . . . . 3.1. The use of Pad4 approximants . . . . . 3.2. The use of rational interpolants . . . . 3.3. Predictor-corrector methods . . . . . 3.4. Numerical results . . . . . . . . . . . 3.5. Systems of first order ordinary differential
. . . . .
. 238 . 241 . 243
. . . . .
. . . . .
. . . . .
. . . . .
. . . . . . . . . . . . . . . . . . . . equations . . . .
238
244
. 247
.
54 Numerical integration . . . . . . . . . . . . . . . . . . . . . 250 4.1. Methods using Pad4 approximants . . . . . . . . . . . . . 251 4.2. Methods using rational interpolants . . . . . . . . . . . . 252 4.3. Methods without the evaluation of derivatives . . . . . . . . 253 4.4. Numerical results for singular integrands . . . . . . . . . . 254
55 . Partial differential equations
. . . . . . . . . . . . . . . . . . 257
§6.Integral equations . . . . . . . . . . . . . . . . . . . . . . . 260 0.1. Kernels of finite rank . . . . . . . . . . . . . . . . . . 260 6.2. Completely continuous kernels . . . . . . . . . . . . . . 262
Problems
. . . . . . . . . . . . . . . . . . . . . . . . . . .
265
Remarks
. . . . . . . . . . . . . . . . . . . . . . . . . . .
267
. . . . . . . . . . . . . . . . . . . . . . . . . .
268
References
196
“It i s my hope that by demonstratdng the e m e with which the various transformations may be effected, their field of application might be widened, and deeper insight thereby obtained into the problems for whoae aoh6tion the tran8formatdons have been uaed.
P . WYNN (1956).
- “On a
device for computing the e m ( & )
tran~f~rmation~’
Z V . l . l . The univariate €-algorithm
197
The approximations introduced in the previous chapters will now be used to develop techniques for the solution of various mathematical problems: convergence acceleration, numerical integration, the solution of one or more simultaneous nonlinear equations, the solution of initial value problems, boundary value problems, partial differential equations, integral equations, etc. Since these techniques are based on nonlinear approximations they shall be nonlinear thernselves. We shall discuss advantages and disadvantages in each of the sections separately and illustrate their use by means of numerical examples.
§I. Convergence acceleration. 1.1. The univariate r-algorithm. Consider a sequence { a ; } ; , ~ of real or complex numbers with lim ai = A
i-
00
Since we are interested in the limiting value A of the sequence we shall try to construct a sequence (b;),.,N that converges faster to A , or
We shall describe here some nonlinear techniques that can be used for the construction of { b;},N. Consider the univariate power series
i-t
with Vai = a; - a;-1. Then clearly for the partial sums
we have
Fk(1)= ak
k = 0, 1,2, . . .
If we approximate f(z) by r ; , ; ( z ) , the Pad4 approximant of order (i,i) for f , then we can put b; = r;,;(l)
i = 0 , 1 , 2 , . ..
198
I V . l . l . The univariate €-algorithm
For the computation of b; the c-algorithm can be used:
Then
The convergence of the sequence {biIiEm depends very much on the given sequence { a , } ; , N . Of course the convergence properties of { b ; } ; , ~are th e same as those of the diagonal Pad6 approximants evaluated a t z = 1 and for this we refer to section 4 of chapter 11. In some special cases it is possible to prove acceleration of the convergence of ( a ; } i E N . A sequence {a;}iEN is called totally monotone if
Ak u ; > O
d,k=0,1,2
,...
where A k a , = Ak-' ai+i - Ak-' ai and Ao a; = a;. In other words, {ai}iEN is totally monotone if
a0 2 a1 2 a2 2 ... 2 0 Aao 5 Aal 5 h a 2 5 . .. 5 0 A2ao 2 A 2 a l 2 A2a2 2 . . . 2 0 and so on.
IV.1.1. The uniuariate 6-algorithm
199
Note t ha t every totally monotone sequence actually converges. The sequences (Xi}i,-~ for 0 <_ X 5 1 and {l/(i l ) } ; , ] ~are for instance totally monotone sequences. The close link with the theory of Stieltjes series becomes clear in the following theorem [4 p.811.
+
Theorem 4.1. The sequence { o ~ } ~ €isN totally monotone if and only if there exists a realvalued, bounded and nondecreasing function g on [0, 11 taking on infinitely many different values such that 1
a; = I 0
t' d g ( t )
i = 0,1,2,. ..
A sequence { a ; } ; , ~is called totally oscillating if the sequence
ai}iE]N
is totally monotone. One can see th at every convergent totally oscillating sequence converges to zero. For these sequences the following results can be proved [4 pp. 83-85]. Theorem 4.2.
If the 6-algorithm is applied to the sequence { U ; } ~ ~ ] N with 1imi-m a i = A and if there exist constants a and @ such that { a a ; + @ } i E ~is a totally monotone sequence, then lim
cfi
=A
lim
($2
k >_ 0 and fixed
=A
e 2 0 and k e d
f-00
k+m
If also lirn
i-oo
then
ai+l
-A # 1
-_-
ai-A
200
IV.1.1. The univariote t-algorithm
Theorem 4.3.
If the €-algorithm is applied t o the sequence { a ; } i E mwith limi-a; = A = 0 and if there exist constants cr and B such that (&a; + @ ) i E ~is a totally oscil-
lating sequence, then
We illustrate these theorems with some numerical results. Consider
This is a totally monotone sequence. In table 4.1. we have listed the values for n = 0 , . . ., 4 and m = 0, . . ., 10. So the values b i = can be found on the main diagonal. Clearly the convergence of the sequence { a i } i e ~is accelerated.
EE)
All these computations were performed using double precision arithmetic (56digit binary representation). To illustrate the influence of data perturbations and ( m - 4 obtained using single rounding errors we give in table 4.2. the same values E~~ precision input and single precision arithmetic (24-digit binary representation). In the tables 4.3. and 4.4. the respective results are given for the totally oscillating sequence
Here the convergence of the sequence {b;}+-N is much faster th a n th a t of the sequence {ai}i,-N. These examples illustrate that the influence of data perturbations and rounding errors can be important. When is nearly equal to one can loose a lot of significant digits. In [62] i t is shown how this misfortune may be overcome.
tP1)
Table 4.1.
56-digit binary reprea ent otio n %
m\n
S 0
1
2
3
4
0
0.10000000D+01
0.86666667D+00
0.52173913D+OO
0.43636364D-bOO
0.37874803D-I-00
1
0.50000000D+00
0.25000000D+00
0.17847058D+00
0.13801453DfOO
0.113Q8283D-tOO
2
0.33333333D-tOO
0.16668667D+00
O.llllllllD~00
0.85271318D-01
0.69561659D-01
3
0.25000000D+00
0.12500000D+00
0.83333333D-01
0.62500000D-01
0.50607287D-01
4
0.20000000D~00
0.10000000D+00
0.66666667D-01
0.50000000D-01
0.40000000D-01
5
0.16666667D+00
0.83333333D-01
0.55555556D-01
0.41666667D-01
0.33333333D-01
6
0.14285714D-kOO
0.71428571D-01
0.47610048D-01
0.35714286D-01
0.2857142QD-01
7
0.12500000D+00
0.62500000D-01
0.41666667D-01
0.31250000D-01
0.25000000D-01
8
O.lllllIllD~OO
0.55555556D--01
0.37037037D-01
0.27777778D-01
0.22222222D-01
9
O.lOOOOOOOD+OO
0.50000000D-01
0.33333333D-01
0.25000000D-01
0.20000000D-01
10
0.9090Q091D-01
0.45454545D-01
0.30303030D-01
0.22727273D-01
0.18181818D-01
k !Y
% F:
2. c P
3.
R
cs
hl
0
N
Table 4.2.
24 -digit binary repreaent ation
m\.
0
1
2
3
4
c..
0
0.10000000E+01
0.86666680E+00
0.52173018E+00
0.436363703+00
0.378748OOE+OO
k
1
0.500000003~00
0.250000083+00
0.17647058E+00
O.l3801455E+OO
0.113082853~00
Y
2
0.33333334E-kOO
0.16666664E+00
0.11111123E+00
0.852712323-01
0.69561742E-01
E
3
0.25000000E+00
0.12500004E+00
0.83333202E-01
0.62500395E-01
0.50605543E-01
4
0.20000000E+00
0.10000000E+00
0.66666812E-01
0.49907089E-01
0.40013745E-01
5
0.16666667E+00
0.83333351E-01
0.55554848E-01
0.416750723-01
0.332807823-01
6
0.14285715E+00
0.71428463E-01
0.47620941E:-01
0.356953373-01
0.28692538E-01
7
0.12500000E+00
0.62500104E-01
0.41664604E:-01
0.312740473-01
0.247051133-01 0.22477847E-01
8
O.llllllllE~00
0.55555556E-01
0.370379803-01
0.27756441E-01
0
0.10000000E+00
0.5000001OE-01
0.33332549E-01
0.2502411OE-01
O.l0564494E:-Ol
10
0.90909004E-01
0.45454517E-01
0.30304417E-01
0.226919743-01
0.18706985E-01
z 2. C
e
2.
f: m
m
it
rg 0 7
3
Table 4.3.
56-digit binary repreaentatran Ls,
m\.
5
0
0
0.10000000D+01
1
-0.50000000D~00
2
0.33333333D.fOO
3
4 5 6
7 8
-0.25000000D+00
1 0.40000000D~00 0.35714286D-01 -0.98039216D-02 0.40322581D-02
0.20000000D~00 -0.20408163D-02 -0.16666667D+00 0.14285714D-kOO -0.12500000D+00 0.111111llD+OO
0
-0.100OOOOOD+00
10
0.~0009001D-01
0.1 17370881)-02 -0.73637703D-03 0.40212508D-03 -0.34506556D-03 0.25125628D-03 -0.18860807D-03
2 0.25531015D-kOO -0.2 18078 10D-01 0.1 1454754D-02 -0.25960540D-03 0.83229297D--04 -0.330401 llD-04 0.15170805D-04 -0.77490546D-05 0.42861980D-05 -0.252506 123)-05 0.1 5651583D-05
3 0.18604651 D+OO 0.11833001D-01 -0.38107008D-03 0.34806706D-04 -0.73003358D--05 0.20314308D-05 -0.688Q106OD-06 0.27006940D-06 -0.11826320D-06 0.58510073D-07 -0.289770043)-07
4
0.14600870D+00 -0.96160448D-02 0.15341245D-03 -0.89086733D-05 0.10459978D-05 -0.20044207D-06 0.53279358D-07 -0.16203117D-07
; Y
z C
2. C P -.
l
a
cs m
a
cq
*
a 0
0.5656006OD-08 -0.220288523)-08 0.0378304OD-00
ts 0 W
01
0
Table 4.4.
b b
m+l Rd-digit binary representation
m n 0
0 0.10000000E~01
1
-0.500OOOOOE~00
2
0.333333348+00
3
4 5
-0.25000000E+00 0.20000000E+00 -0.16666667E+00
6
0.142857153+00
7
-0.12500000E~00
8
9 10
O.llllllllE~00 -0.10000000E+00 O.QO909094E-01
1
0.400000043+00 0.357142843-01 -0.880381823-02 0.40322733E-02 -0.20408179E-02 0.11736Q87E-02 -0.738376743-03 0.48213349E-03 -0.345061333-03 0.25125383E-03 -0.18860842E-03
2 0.25531918E-kOO -0.218878OOE-01 0.1145486OE-02 -0.25958874E-03 0.832260823-04 -0.3305387QE-04 0.15181173E-04 -0.77439190E-05 0.428871413-05 -0.25255961E-05 0,156517563-05
4
3 0.18604654E3+00 0.118331053-01 -0.38106500E-03 0.34898013E-04 -0.73038805E-05 0.20300561E-05 -0.68641026E-06 0.273580183-08 -0.11687172E-08 0.56633620E-07 -0.28493082E-07
0.146908823+00 -0.9816957OE-02 0.15340716E-03 -0.881250783-05 0.10442981E-05 -0.20828127E-06 0.56101026E-07 -0.138258973-07 0.860953913-08 -0.17431754E--08 0.183647753-08
* Y
k k
Y
z E
a.
a 3. C
h R
m
$-
-l
3
I V . 1.2. The qd- algorithm
205
Various generalizations t o accelerate the convergence of a sequence of vectors, matrices, and so on exist. We refer to [4] and [Sl]. The convergence of a multidimensional table of real or complex numbers is treated in section 1.5. 1.8. The qd-algorithm.
To accelerate the convergence of { ai}iEN we construct again
i- 1
We have seen in section 3 of chapter 11 th at it is possible t o construct a corresponding continued fraction I
I
I
where the coefficients d i can be calculated by means of the qd-algorithm. If we calculate the i t h convergent Ci(z) of this continued fraction, we can now put 6i = Gi(1)
i = 0 , 1 , 2 , .. .
The numbers 6; can be computed using a continued fraction algorithm from section 5 of chapter I. Remark that theoretically the numbers b 2 ; here are the same as the 6; computed by means of the €-algorithm. In practice there are however differences due to rounding errors. We will illustrate this numerically. The tables 4.5. and 4.6. illustrate the effect of perturbations when the qd-table is computed: the results in table 4.5. are obtained using double precision arithmetic while those in table 4.6. are obtaincd using single precision arithmetic. Input was
So f(z)= ez and limi-m a; = e = 2.718281828.. . From the tables 4.5. and 4.6. the values b; were computed using the three-term recurrence relations (1.3.) respectively in double and single precision. The results can be found in the respective tables 4.7. and 4.8., where they can also be compared with the values of the c-algorithm.
Table 4.5.
hl
0 Q,
k=l,2,
92
-0.16668667D-kOO 0.33333333D+00
-0.23800524D-01 0.142857 14D+00
-0.17857143D-01 0.12500000D+00
0.97222222D-01
-0.1388888OD-01 0.1 11111 11D+OO
-0.90900091D-02 0.90000001D-01
0.727272731)-01
0.18181818D-01
0.81818182D-01
-0.27272727D-01
-0.45454545D-01
0.54545455D-01
-0.36363636D-01
0.63636364D-01
(4
e5
0.55555556D-01
-0.44444444D-01
0.66666667D-01
-0.33333333D-01
0.77777778D-01
-0.22222222D-01
0.8888888OD-01
-0.11111111D-01 0.10000000D+00
-0.27777778D-01
(kf
95
-0.55555556D-01
0.60444444D-01
-0.41666667D-01
0.83333333D-01
(a)
e4
0.71428571D-01
-0.53571429D-01
0.80285714D-01
-0.35714286D-01
0.107 14286D-kOO
(k)
94
-0.7142857lD-01
0.05238005D-01
-0.47610048D-01
0.110047821)+00
(k) e3
0.10000000D+00
-0.66666667D-01
0.13333333D-kOO
-0.33333333D-01 0.16866867D+00
(k)
93
-0.10000000D+00
0.15000000D~00
-0.50000000D-01 0.20000000D+00
(k)
e2
O.I6666667D+OO
-0.83333333D--01 0.25000000D+00
56-digit binary reprerent ation
(k)
0.50000000D~00
...
Y t r
.4
n
Table 4.6. gd - table f o r gik’ =
Vak+1
v ak
_ I _
k
with ak =
C -j!1
k = 1,2, ...
j-0
4 - d i g i t binary representation
0.333333343+00
-0.833333433-01 0.25000000E~00
0.16666666E+00
0.14285713E-l-00
-0.17857134E-01 0.12500000E~00
-0.13888896E-01 0.111 llllOE+OO
0.888887723-01
-0.11111103E--01 0.100000003+00
-0.454895578-01
0.545119843-01
-0.36346123E-01
0.836478143-01
-0.272756813-01
0.72725780L?3-01
-0.1 81818443-01
0.818182303-01
-0.0090007OE-02 0.80809084E-01
-0.222224373-01
(k)
e5
0.555540553-01
-0.444474973-01
0.666824538-01
-0,333314843-01
0.77779025E-01
(4
95
-0.55559866B--01
0.694431893:-01
-0.416667613-01
0.83332881E-01
-0.277776273-01
0.972223 13E-01
(k)
e4
0.714307433-01
-0.535723123-01
0.89285374E-01
-0.357142693-01
0.107142813$00
(i)
94
-0.714278823-01
0.052385073-01
-0.47619178E-01
O.l1804755E+OO
-0.238085223-01
(4
e3
0.99Q89800E-01
-0.666665883-01
0.13333338E+00
-0.333333483-01
93
-0.10000007E+00
O.l4999998E+OO
-0.49999989E-01 0.20000000E+00
(k)
0.16666670E+00
IV.1.2. The qd-algorithm
208
Table 4.7.
lim a; = 2.71828182845904..
i-+co
j-0
56-digit binary repreaentation (0)
i
€2 i
1.000000000000000D+00 1.000000000000000D+OO 3.000000000000000D+OO 3.000000000000000Df00 2.7142857142857 14D+00 2.714285714285714D+OO 2.7 18309859154930D+00 2.718309859 154930D+00 2.718281718281718D+OO 2.718281718281718D+OO 2.7 18281828735696D+00 2.718281828735696D+OO 2.7182818284585631)+00 2.7 18281828458564D+00 2.718281828459046D+OO 2.718281828459046D+OO 2.7 18281828459045D+00 2.718281828459045D+OO Table 4.8. 1
a i = j -C 0,
lim ai = 2.718281828.. .
i-
00
24-digit binary represent ation 1
0 1
2 3 4 5
I 1.0000000E+00 3.0000000E+00 2.7142858E+00 2.7183099E+00 2.7 1828183+00 2.7182815Ei-00
1.0000000E+00 3.0000000E+00 2.7 I42858E+OO 2.71830993+00 2.71828183+00 2.71828223+00
.
IV.1.3. The dgorithm of Bulirach-Stoer
209
1.3. The algorithm of Bulirsch-Stoer. Let
(2,)
i
e be ~ a convergent sequence of points with
limi-+w z i = z When using extrapolation techniques to accelerate the convergence of we compute a sequence {bi};,N with
{U;},.~N
b; = limz-.z gi(z) where g i ( z ) is determined by the interpolation conditions gi(zi) = a j
j
0 , . . ., 9: with i
2
0
The point z is called the extrapolation point. Often polynomials are used for the interpolating functions 8 ; . We consider here the case of rational interpolation where we shall use Stoer’s recursive scheme given in section 3 of chapter 111. Using the notations of chapter I11 and writing
we first generalize the formulas (3.12.):
and
with pO,o(s) ( i ) = a, and qo,o(z) (i) = 1. We can easily see taking into account the interpolation properties of r(J+l)
m,n- 1
rk!n and
that
(4.1a .)
IV.1.3. The algorithm of Bulirech-Stoer
210
and analogously that
(4.1b .)
If aiT;ll
# 0, then we can write
or
Analogously
If we use the interpolating functions
IV.f.3.The algorithm of Bulirach-Stoer
and write 8 (; )' = r@,(z) for rn
with
= t-tb(z)= a j and
211
+ n = &,then
8-1 ( j1
= 0 for j = 0 , 1 , 2 , . . , and
Let ua compare t h e algorithm of Bulirsch-Stoer with the €-algorithm, the qdalgorithm and the method of Neville-extrapolation based on polynomial interpolation (here bi = rri(z)). Input is the sequence {ai);,N = { ~ / r n ) ; ~ N
with lim;+a a; = 0. For the methods based on interpolation we use
z=
lim x i = O
i- oa
IV.i.3. The algorithm of Bulirsch-Stoer
212
The results are listed in table 4.9. Computations were performed in single precision arithmetic.
Table 4.9. a.=
1
-__
&Ti Bulirsch-
i ~
qd
Neville
Stoer
epsilon
Qi _____--_-- b i
s(0)
€ i(0)
C;(1)
0
1.00000000
1.00000000
1.00000000
1.00000000
1.00000000
2
0.5773502e
0.26964909
0.24118084
0.47414386
0.47414389
4
0.4 472 13 59
0.16024503
0.11 8 2 0 6 8 3
0.3095461 1
0.3095461 1
6
0.3779e450
0.11422185
0.07439047
0.22860453
0.22960511
8
0.33333334
3.08880618
0.06251375
0.18284784
0.18285738
10
0.30151 1 3 5
0.07302472
0.05196425
0.15728275
0,15730451
12
0.27735010
0.0617 1018
0.04684320
0.14 1 I5 857
0.14100026
14
0.25819880
0.11071037
0.04458576
0.12546575
0.12571160
0.04921301
0.10151753
0.10491328
0.34646210
16
0.24253564
I8
0.22941573
-2.86674523
0.04912448
0.10688034
0.10763902
20
0.2 1821788
-89.92565918
0.05485482
0.09624152
0.10446023
I V . l . 4 . The p-algorithm 1.4. The
2 13
p-algorithm.
Another technique is based on the construction of interpolating continued fractions of the form
If we take for g;(z) the ith convergent
and take z = 00 as extrapolation point, then we can put for i even:
6; = lim 30’2
g i ( z ) = do
+ d:! + . . . + di
Now the sequence b i can be calculated using reciprocal differences. Hence we can set up the following scheme:
k =0,1,2,. .. Zk-1. k = 1,2,3,... tk,l = a k - ak-1 Xk - Z k - j __ tk,j = tk-l,j-2 + tk,j-l - tk-1,j-1 tk,O
= ak
Zk
k=j,j+l
where in fact, in the notation of chapter 111, t k $ 3' = p [ Z k - j ,
. . .,zk]
Then clearly for a even
.,
bi = ~ [ z o , . .z;]
The values
tk,j
= ti,i
are ordered as in the following table
,...
and j = 2 , 3
,...
I V . l . 4 . The p-algorithm
214
to ,o t l ,o
tl,l
t2 ,o
t2,l
t2,2
t 3 ,O
t3,l
t3,2
with t k , o = a k and t k , - i = 0 for k = 0 , 1 , 2 , . . .. if we use { z ; } ~ ~as N interpolation points for the algorithm of Bulirsch and Stoer and {Z;};~IN with
as interpolation points for the p-algorithm, then for even 4:
(see problem (3)). Remark t ha t the computation of t i , i takes a smaller effort than th a t of 8y). For more properties of the p-algorithm we refer to [4]. We illustrate the influence of the choice of the {z~},.~N with some numerical examples. In the tables 4.10. and 4.11. several { b i } i E m were constructed, respectively for {ai);,N
={ l / m } i E N
and {a;)ieN = {I/('
+l ) } a ~
All computations were performed in single precision arithmetic. In each table the first sequence {b;}i,m is constructed with xi = &,the second with z i = i, the third with Z; = i2 and the fourth with xi = e i . The choice of the 2; greatly influences the convergence behaviour of the resulting sequence {b;};,N. We see that we do not necessarily obtain better results when the 2 ; converge faster to infinity. More information on this matter can be found in [8].
215
I V . l . 4 . The p-algorithm
Table 4.10.
___ i
-t ,. b; = t i , ; 6. b; 2 ti,; 1 - *,I 2; = i 2 z’- -2; = i a; I -_ _ ~ -
6
b-=t., x; = e' 1
1,t
0
1.00000000
1.00000000
1.00000000
1.oooooooo
1.00000000
2
0.57735026
7.07813832
0.24 11 8 0 8 8
0.50412308
0.48505624
4
0.44721358
0.02822608
0.11828871
0.30880488
0.30424682
6
0.37786450
-0.00126043
0.07431557
0.22015348
0.34178661
8
0.33333334
-0.00124686
0.06307473
0.17057084
0.30704348
10
0.30151135
-0.004163Q2
0.05177462
0.1 4054570
0.28143203
12
0.27735010
0.04758168
0.12608825
0.26141584
14
0.25818889
0.04435266
0.1 1 4 5 1 1 1 3
0.24517243
16
0.24253564
0.00017822
0.04088086
0.10530381
0.23163427
18
0.2284 1 5 7 3
0.00026065
0.05084635
0.00724047
0.22012024
0.21821788
0.00000440
0.03507820
0.08777263
0.21018887
20
0.00010506 -0.00047068
Table 4.11. a; =
i __
b*.
ai
0
1.00000000
2
0.33333334
4
0.20000000
6
0.14285716
8
0 . 1 1111111
10
0.08080008
=z
1 i+l
t.. 1,'
xi = -
b; = t %. *. bi = t i , ; .2 2; = i -~ --xi = 9,
1.00000000
b . = t *. , i
xi = e'
1.00000000
1.00000000
0.00000006
0.25000003
0.2401564 1
0.00000162
0.00000003
0.08080016
0.1543827 8
0.00000038
0.00000002
0.04545450
0.31652598
-0 .OOOOOO 15
0.02703313
0.08417003
0.00000003
0.01783462
0.07815881
-2.41421294
-0.00000042 0.00000018
1.00000000
IV.1.5. The multivariate e-algorithm
216
In [49] a study was made of some nonlinear convergence accelerating methods. The +algorithm was found to be the most effective one in many cases. Therefore we now generalize the t-algorithm for the convergence acceleration of multidimensional tables. 1.5. The multivariate t-algorithm. can he considered as a table with single entry. For its convergence acceleration we constructed the univariate function
A sequence {ai},.,N
c 00
f(z) =
VUiX'
i=O
with ai = 0 for i
< 0 and calculated
diagonal Pad6 approximants evaluated a t
x = 1. Let us now first consider a table (aj,k
}j,ken
with double entry and with A = limj,k.+OOaj,k = limj-w(limk+OO ajk) = limk-+w(limj+m ajk). To accelerate its convergence we introduce
We will again construct a sequence ( b i } i E ~ that will converge faster t o A in some cases [ll].To this end we now calculate bivariate Pad6 approximants for f ( z , y) and evaluate them at ( 5 , y) = ( 1 , l ) . If we denote by Ut =
j+k==t
Qjk
= at,O
+ at-1,l + . . + al,L-1 + a0,L
and start the c-algorithm with th e partial Bums of f( 1 , l )
IV.1.5. The multivariate 6-algorithm
approximants ri,i(z,y) evaluated a t
then the bivariate diagonal l’adi: (0) (e, y) = (1,l) are given by t2i . Hence we put bi
217
(0) = e2i
Let US generalize the idea for a table with multiple entry define
{ ~ j ~ . , , j ~ } ~ ~ We , , . . , ~ , ~ ~
W
f ( z l , .. .,zk) =
jl
, .. . ,j,=O
Vajl.,.j, Zil.
..zc
with
It is easy to prove that
Again multivariate Pad6 approxinlants for f ( z l , . . ., zk) can be calculated and evaluated a t (z,, . . .,zk)= ( 1 , .. ., 1) via the t-algorithm. Sirice
where now
the
6f)
for f(z1,. . .,Zk) are given by
IV.1.5. The multivariate €-algorithm
218
We illustrate this type of convergence acceleration with a numerical example. Suppose one wants to calculate the integral of a function G ( z ~ , ..,. Zk) on a bounded closed domain D of Rk. Let D = [0,1] x . . . x [0, I] for the sake of simplicity. The table { a i l . . . i k } j ~ , , , , , j can k ~ ~ be obtained for instance by subdividing the interval [O,1] in the lthdirection (t! = 1 , . . ., k) into 2jl intervals of equal length h~ = 2-Je(jt = 0, 1 , 2 , . . .). Using the midpoint-rule one can then substitute approximations
to calculate the ail...j, . We restrict oursclves again t o the case of two variables. With hl = 2-j and hp = 2-k we get
The column c0(4( l = 0, 1, . . .) in the 6-table given by
was also used by Gens [21] to start the €-algorithm for the approximate calculation of multidimensional integrals by means of extrapolation methods. He preferred this method to six other methods because of its simplicity and general use of fewer integrand evaluations. For
we have
[1
1
lim
j,k-cce
ajk =
1
dzdy = 2 In 2 = 1.386294361119891...
In table 4.12. one can find the a i , slowly converging to the exact value of the integral because of the singularity of the integrand in the origin. The 6; = r,,;(l, 1) converge much faster. For the calculation of the b i we need the c~ ( l = 0 , . . ., 2t).
IV.1.5. The multivariate €-algorithm
Hence b; should be compared with an information, i.e. j -tk = 2i.
ajk
which uses the same amount of
Table 4.12.
d
aii
1 2 3 4
1.166666606667 1.269017019048 1.325743700744 1.355532404415
219
bi = t2; (0)
1.330294906166 1.396395820203 1.386872037696 1.386308917778
220
IV.2.1. Iterative methods brrsed on P o d i approzimation
$2. Nonlinear equations.
Suppose we want to find a root xu of the nonlinear equation
f (4 = 0
Here the function f may be real- or complex-valued. If f is now replaced by a local approximation then a zero of th at local approximation could be considered as an approximation for x*. Methods based on this reasoning are called direct. One could also consider the inverse function g of f in a neighbourhood of the origin, if it exists, and replace g by a local approximation. Then a n evaluation of this local approximation in 0 could be considered as an approximation for x* since g(0)
.=z *
Methods using this technique are called inverse. 2.1. Iterative methods baaed on Pad6 approzimation.
Let
2;
be an approximation for the root xu of ri(x)
and let
Pi
= -(x) cli
be the Pad6 approximant of order ( m ,n) for f in z i . Then the next approximation x;+1 is calculated such that Pi(Zi+l)
=0
In case p ; ( z ) is linear ( m = 1) the value zi+l is uniquely determined. It is clear that this is t o be preferred for the sake of simplicity. A well-known method obtained in this way is Newton’s method (rn = 1, n = 0) which can be derived as follows. The Taylor series expansion for f(z) at zi is given by f(5) = f ( 5 ; ) f f ’ ( Z i ) ( Z
- 5 i ) + f‘ I ( x i ) ( 5 - X i ) 2 2
Hence the Pad6 approximant of order (1,O) for f a t
xi
+. . *
(4.2.)
equals
and we obtain (4.3.)
IV.2.1. Iterative methods boaed on Pad4 approzimation
221
Another famous method is Halley’s method based on the use of the Pad4 a p proximant of order (1,l) for f at X i : (4.4.)
Since the iterative procedures (4.3.) and (4.4.) only use information in the point z; to calculate the next iteration point z;+1 they are called one-point. The order of convergence of these iterative methods is given in the next theorem.
Theorem 4.4. If { z ~ } ~converges ~ N t o a simple root z* of f and if r i ( z ) is normal for every i then
For the proof we refer to [52]. Similar results were given in [IS], [43] and 1441. So the order is a t least rn n 1 and it depends only on the sum of m and ra. Consequently the order is not changed when Pad6 approximants lying on a n ascending diagonal in the Pad4 table are used. Iterative methods resulting from the use of (m,n) Pad6 approximants with n > 0 can be interesting because the asymptotic error constant C* may be smaller than when n = 0 [43]. The iterative procedures (4.3) and (4.4) can also be derived as inverse methods (see problem (5)). Let us apply Newton’s and Halley’s method t o the solution of
+ +
The root Z* = 0. We use xo = 0.09 as initial point. The next iteration steps can be found in table 4.13. Computations were performed in double precision accuracy (56-digit binary arithmetic).
222
IV.8.1. Iterative methods baaed on P a d 4 approzimation
Table 4.13.
i
Newton
Halley
Zi
2;
0
O.QOOOOOO0D-01
1
0.80997588D-01
2
0.6559965QD-01
3
0.43022008D-01
4
~.18500311D-01
5
0.34212128D-02
6
0.11703452D-03
7
0.13697026D -06
8
0.18760851D-12
9
0.35197080D-24
10
0.00000000D+00
0.90000000D -0 1 -0.40647645D-02 0.35821646D-06 -0.24514986D-18 0.00000000D+00
It is obvious that a method based on the use of (m,n) Pad4 approximants for f with n > 0 gives better results here: the function f has a singularity at z = 0.1. Observe that in the Newton-iteration zg is a good initial point in the sense that from there on quadratic convergence is guaranteed: 1zi+1
- 2.1 = 1zi+l I N c ' J z ~- z * I 2 for i 2
o
with C’ = 10. For the Halley iteration we clearly have cubic convergence from X I on. The formulas (4.3.) and (4.4.)can also be generalized for the solution of a system of nonlinear equations
which we shall write as
IV.d.1. Iterative methods based on Pad6 approzimation
223
Newton's method can then be expressed as [45]
where F’(xr),. . . ,xr)) is the Jacobiau matrix of first partial derivatives evaluated at (zy),. .. ,zt))with
Let us now introduce the abbreviations Fi = F ( z f ) ,. . . ,xk(i))
F , ! = F ( z (4 , ,... ,zk (4) To generalize Halley's method we first rewrite (4.4.) as
Then for the solution of a system of equations it becomes [I41
224
IV.2.1. Iterative method8 based o n P a d 6 appton'mdion
where the division and the square are performed componentwise and F " ( x 1 , ... ,x k ) is the hypermatrix of second partial derivatives given by
~~
a 2f k a2f k ... a2f k axkaxs ... axlaxk
dxlax2
.
.
I
az:
which we have t o multiply twice with the vector --F,!-'Fi. This multiplication is performed as follows. The hypermatrix F " ( x 1 , . . . , Z k ) is a row of k matrices, each le x k. If we use the usual matrix-vector nlultiplication for each element in the row we obtain a2fi
c
a2fk
c
k
asl ax;Yi ... i= 1 axlax;
Yi
...
k
i= 1
jin) vk
In [14] is proved that the iterative procedure (4.6.) actually results from the use of multivariate Pad6 approximants of order ( 1 , l ) for the inverse operator of ~ ( x l. ., . , x k ) at ( x v ) , . . . , x g ) ).
IV.8.1. Iterative method8 baaed on Pod6 approzimation
225
To illustrate the use of the formulas (4.5.) and (4.6.) we shall now solve the nonlinear system
{ (-
f 1 (‘J
Y) = -=+I/ = e-=-y
f&,y)
- 0.1
=0
- 0.1 = 0
which has a simple root at 12.0.1)) = (2.302585092994O46.. 0.
.
As initial point we take ( d 0 ) , g ( ’ ) ) = (4.3,2.0). In table 4.14.one finds the consecutive iteration steps of Newton’s and Halley’s method. Again Halley’s method behaves much better than the polynomial method of Newton. Here the inverse operator G of the system of equations F has a singularity near to the origin and this singularity causes trouble if we get close to it. For
we can write -0.5(ln(0.1 0.5(ln(0.1
+ u) + ln(O.l + v)) + u) - ln(O.l + v))
With (do), y(O)) = (4.3,2.0) the value do)= f2(z(O),y(O)) is close t o -0.1 which is close t o the singularity of G. For the computation of th e Pad6 approximants involved in all these methods the €-algorithm can be used. Another iterative procedure for the solution of a system of nonlinear equations based on the e-algorithm but without the evaluation of derivatives can be found in [5]. Since it does not result from approximating the multivariate nonlinear problem by a multivariate rational function, we do not discuss it here.
a
J9 0
Ne ton
0.43000000D+01
H Iev
,A4 0.20000000D+01
N N
0.43000000D+01
0.20000000D+01 0.5705728dD-kOO
1
-0.22427305D+02
-0.24720886D+02
0.287Q8400D+01
2
-0.21927303D+02
-0.24228888D-k.02
0.22495475D+Ol
-0.52625816D-01
3
-0.21427303D+02
-0.23720888D-kO2
0.23018229D+Ol
-0.44Q01947D-02
4
-0.20827303D+02
-0.23228888D+02
0.23025841D+01
-0.57154737D-05
5
-0.20427303D+02
-0.2272Q888D-kO2
0.23025851D-l-01
-0.Q7680305D--11
6
-0.1 QQ27303DS-02
-0.22228888D-kO2
0.23025851D+Ol
-0.17438130D--16
7
-0.18427303D+02
-0.21720888D-bO2
8
-0.18Q27303D+02
-0.2122Q888D-tO2
Q
0.18427303D+02
-0.20728888D)+02
10
-0.17827303D)+02
-0.20229888D-kO2
11
-0.17427303D+02
-0.19720888D-kO2
12
-0.16927303D+02
-0.10228888D+02
13
-0.16427303D+02
-0.18728888D+02
14
-0.15Q27303D+02
-0.18228888D-I-02
15
-0.15427303D+02
-0.17729888D+02
16
-0.14927303D+02
-0.17229888D+02
17
-0.14427303D+02
-0.16729888D+02
18
-0.13027303D+02
-0.1622Q888D-tO2
18
-0.13427303D+O2
-O.l572Q888D+02
20
-0,12!327303D+02
-0.15228888D-tO2
rn
y(4
3 e
m
J
?:
;r a
R e
Iv.8.B. Iterative methods based on rational interpolation
227
3.2.Iterative methods baa e d o n rut i o n a1 int erpol at i o n. Let
with pi and 9, respectively of degree m and n, be such that in an approximation z i for the root z * of f
fr)(zi-j)
with m that
= f(')(Zi-j)
I!
= 0 , . . .,8j - 1
(4.7.)
-+ n + 1 = C'tz08 ~ Then . the next iteration step z;+] is computed such Pi(zi4-1)
=0
For the calculation of zi+l we now use information in more than one previous point. Hence such methods are called multipoint. Their order of convergence can be calculated as follows. Theorem 4.5.
If
(Z;};~Nconverges to a simple root z* of f and f(n+n+l)(z) continuous in a neighbourhood of z* with
with n
> 0 is
where f ( k ) ( z *= ) 0 if k < 0, then the order of the iterative method based on the use of r,(z) satisfying (4.7.) is the unique positive root of the polynomial .i+1
- 8ozj- B1zj-l
The proof can be found in 1521.
- .. .- 8j = 0
228
IV.2.2.Iterative methods baaed on rational interpolation
If we restrict ourselves to the case 8[
=8
l = 0 , . , ., j
then it is interesting t o note that t h e unique positive root of
!=O
increases with j but is bounded above by 8 + 1 [53 pp. 46-52]. As a conclusion we may say t h a t the use of large j is not recommendable. We give some examples. Take m = 1, n = 1, 8 = I and j = 2. Then z;+1 is given by
The order of this method is 1.84, which is already very close to 8 + 1 = 2. Take m = 1, n = 1, 80 = 2, 81 = 1 and j = 1 . Then x;+1 is given by
The order of this procedure is 2.41. The ease m = I , n = 0,8 = 1 and j = 1 reduces to the secant method with order I .82. Let us again calculate the root of
with initial points close t o the singularity in x = 0.1. The successive iteration steps computed in double precision (56-digit binary arithmetic) are shown in table 4.15.
IV.2.2. Iterative method8 based on rational interpolation
229
Table 4.15.
(4.8.)
(4.0.)
secant method
Xi
zi
Xi
d
0.80000000D-01
0
1 2 3
0.90OOOO0OD-01
0.80000000D-01
0.80000000D-01
0.85000000D-01
O.QO0OOOOOD-01
O.QOOO00OOD-01
-0.15847802D-03 0.15131452D-06 -0.51538368D-13
4
0.20630843D-24
5
0.00000000D+00
- 0.17324153D-03
0.71095917D-01
0 . 4 6 6 2 1186D-10
0.64701421D-01
-0.6139231 1 D-25
0.46636684D-01
0.00000000D+00
0.30206020D-01 0.14080364D-01
6
0.425123451)-02
7
0.59843882D-03
8
0.28 4 30 08 1 D -0 4
0
0.15223576D-06
10
0.38727362D--10
11
0.58956894D-16
12
0.22832450D-25
By means of the multivariate Newton-Pad4 approximants introduced in section 8.2. of chapter 111 the previous formulas can be generalized for the solution of systems of nonlinear equations. We use the same notations as in chapter 111 and as in the previous section. For each of the multivariate functions f,(q,.. .,zk) with j = 1 , . ..,k we choose
D = N = ((0,. . . , O ) , ( l , O , . . .,O ) , ( 0 , 1 , 0 , .. . , O ) , . .., (0,...,o, 1)) H={(2,0,
..., 0 ) , ( 0 , 2 , 0,...10),..., (0,..., 0,2)) 2 INk
c Nk
IV. 2.2. Iterative methods baaed on rational interpolation
230
Here the interpolationset N U H expresses interpolation conditions in the points
zp)). . ., ($',
..
.)
Remark that this set of interpolation points is constructed from only three successive iteration points. The numerator of
with possible coalescence of points, is then given by
NO ,...,O ( 21, . . ,zk) c02,00,. ..(00
coo, ...,00,02
where
Nl,O
,...,0 (21> .
1
* 3
c12,00,...,00
0
zk
. . . NO ,...,0,1(21, . . ... 0 ...
*
coo, ...,00,12
>
zk)
IV.B.8. Iterative methods based o n rational interpolation
231
The values e,,tl,...,,Ltk are multivariate divided differences with possible coalescence of points. Remark that this formula is only valid if the set H provides a systern of linearly independent equations. The next iterationstep (21 ( i + l ), . . ., zk (i+l)) is then constructed such that
p i k ( z y l ) ., .
.J
Z k( i + l ) ) = 0
For k = 1 and without coalescence of points this procedure coincides with the iterative method (4.8.). With k = 2 and without coalescence of points we obtain a bivariate generalixation of (4.8.). Let us use this technique t o solve the sgslem
{
e--l+Y e-Z-9
= 0.1 = 0.1
with initial points (3.2, -0.95), (3.4, -1.15) and (3.3, -1.00). The numerical results computed in double precision (56-digit binary arithmetic) are displayed in table 4.16. The simple root is (2.302585092994046.. . ,O.). In this way we can also derive a discretized Newton method in which the partial derivatives of the Jacobian matrix are approximated by difference quotients
N = H = ( ( 0 , ...,Q ) , ( l , O , ..., 0 ),,..J(Ol...lO,l)}
D
= ( ( 0 , . . .,n))
If we call this matrix of difference quotients AFi, then the next iterate is computed by meana of
As an example we take the same system of equations and the same but fewer initial pointa as above. The consecutive iteration steps computed in double preciaion (56-digit binary arithmetic) can now be found in table 4.17.
232
IV.2.2. Iterative method8 baaed on rational interpolation Table 4.16.
i
~
y")
~-
0.320000OOD+Ol
-0.Q50000OOD-l-00
0.34000000D-k01
-0.11500000Df01
0
0.33000000D-k01
-0.10000000D~01
1
0.25249070D-kOl
-0.22072875D-kOO
2
0.22618832D-tOl
3
0.231278OQD)+Ol
-0.101644QOD-01 -0.51269373D-03
0.41 9 7 l Q 4 4 D - 0 1
4
0.23030978D-kOl
5
0.23025801D-kOl
0.40675854D-05
6
0.28025851Df01
-0.2568682QD-08
7
0.23025851 D+Ol
-0.1277881 6D-13
8
0.23025851Df01
-0.1 1350Q32D-16
Table 4.17.
i 0 1
~~
0.34000000D-k01
-0.11500000D-k01
0.33000000D-b0 1
-O.lO0OOOOOD+Ol
-0.29618530D-kOO
0.21743833Df01 0.20884Q33D-kOl
2
0.32743183D-l-01
3
0 . 2 2 1 1 4 2 1 1 Di-0 1
-0.84011352D+01
4
0.3651533RDt-0 1
-0.72149651D-kOl
5
-0.17Q00083D-kO4
0.208541 11 D i - 0 4
divergence
IV.2.3. Iteratave method3 wing continued fractions
233
The rational method is again giving better results. Now the initial points are
such t h at u = f l ( z , y ) is close to -0.1 which is precisely a singularity of the inverse operator for the considered system of nonlinear equations. For a more stable variant of thc diacretized Newton method we refer to [24]. 2.3. Iterative method8 wing continued fractions.
If ri(z) is t h e rational interpolant of order (m, 1) for $(z), satisfying 1 ri(zw)) = --(&-L))
L = 0,.. . , m + 1
f
then ri(z) can be written in the form rj(z) = do
+ d l ( z - z(i-m-')
+ ldm(z
z(i-m-
-
) + . . . + dm-l(Z - z(i-m-1) 1)
). . .(. - z(i-2J)l
t
+I
2
-
1.. .(z- z ( i - 3 ) )
&-q
&+I
The coefficients d , ( j = 0 , 1 , . . ., m ) are divided differences while dm+l is an inverse difference. The root of ,!;(z) can be considered as an approximation for the root z* of f . SO z(i+l)
-2
(i-1)
- dm+l
This method can be compared with methods based on the use of rational interpolants of order (1, rn) for f(z).
2.4. The qd-algorithm. First we s t a t e an important analytical property of the qd-scheme. To this end we introduce the following nokition. For the function f(s) given by its Taylor series development f(.) = co c 1 2 c 2 22
+
we define the H ankel-det erminant s
+
+
234
IV.Z.4. The gd-algorithm
We now call the series f(z)ultimately k-normal if for every R with 0 5 n 5 k there exists an integer M,, such t h a t for m 2 M,,the determinant Ifm,,,is nonzero.
Theorem 4.8. Let the meromorphic function f(z) be given by its Taylor series development in the disk B(0,r) = {z E C 121 < r ) and let the poles wi of f in B(0, r ) be numbered such that
1
0
< Iw] L
Iwzl
I ... < r
each pole occuring as many times as indicated by its order. If the series f(z) is ultimately k-normal for some integer k > 0, then the qd-scheme associated with this series has the following properties: a) for each n such that 0 < n 5 k and such t h a t Iwn-~l < (tun( < Iw,,+l(, where wg = 0 and, if f has exactly k poles, wk+l = 00, we have Iim m-+m
b) for each n such t h a t 0
1
qLm) = Wn
and such that Iw,(
Iim
m-cc
< Iw,,+lI,
we have
eim) = O
The proof can be found in [29 pp. 612-6131. As a consequence of this theorem the qd-algorithm is an ingenious tool t o determine the zeros of polynomials and entire functions. For if p(z) is a polynomial then the zeros of p ( z ) are the poles of the meromorphic function f(z) = b(z). AIIYq-column corresponding to a simple pole of isolated modulus would tend to the reciprocal value of that pole. It would be flanked by e-columns that tend t o zero. If moreover f ( z ) is rational the last e-column would be zero. Unfortunately the qd-scheme as generated in section 3 of chapter TI is numerically unstable. Rounding errors play an important role due to the divisions by small e-values. However, the rhombus rules (2.6.) and (2.7.) defining the scheme may be rearranged as follows:
IV.2.4. The qd-algorithm
235
In this way quotients are formed only with the quantities q p ) which do not tend t o zero and the qd-scheme is generated row by row. This is called the progressive form of the qd-algorithm. The problem of getting started is solved when we also calculate q y t l ) and elk+’) for negative values of k, i.e. if we calculate the extended qd-scheme as described in section 3 of chapter TI. In [XY pp.641-6541 it is also shown how to deal with the case in which several poles have the same modulus. To conclude this section we illustrate the preceeding theorem with the values qim) and el") for f(z) =
sin x 0.1 - x
which is an ultimately 1-normal function. Clearly with wo = 0,
w1
= 0.1 and
1 lirn qIm) = 10 = Wl
m-w
lim
m-oo
elm)
=0
Table 4.18. 0.10000000000000E+02
e
0 . Q Q S333 33 3 3333 3 E 4-0 1
ey)
=
0 . 1 ~ ~ ~ 6 6 6 ~ 6 6 ~ 6 ~ 7 ~ - n i
o.iooooooooooonnrc-t-o2
el")
=
0.83472454091016E-05
0.10000008347245E~+02
t:
=
-n.s3172454ooioie~-o5
0.100000000000001.;+02
$1
=
-0.108743778QQi25E-08
O.QOOOOQQQQS0l26E+O1
ey)
=
0.10000000000000E+02
e?
O . ~ O ~ O O O O O O O O O O ~ E + ~t:?) ~ 0.10n00000000000E-t02
-0.16686866666687E-01
-
=
0.1Q8743776QQ125E-08 0.27600144302181E-12 -0.276001 44382181 E-12
IV.R.5. The generalized qd-algorithm
238
2.5. The generalized qd-algorithm. The function f(s)can also be given by its Newton series 00
f(z)=
C f[zo, -
7
~ i ] ~ i ( z )
i=O
where the ~ [ z o. ., . , q] are divided differences with possible coalescence of points and where
ll j=l (z - zj-1) i
Bi(s) =
If we define the generalized Hankel-determinants
with Rm,o= 1 then we call the Newton series ultimately k-normal if for every n with 0 5 n 5 k there exists an integer Adn such that for m 1 Mn the determinant Hm,nis non-zero.
Theorem 4.7. Let the sequence of interpolation points ( 5 0 , 5 1 , 5 2 , .. .>be asymptotic to the sequence { ZO, zl,.. . , z,, 20, z1, . . . ,z i , 2 0 , z1, . . . ,z j , . . .} in the sense that
lim
k-oo
sk(j+l)+i =
d = 0 , . . . ,j
z;
Let the function f(s) be meromorphic in
..
~ ( 2 0 , ., Z j , r )
= {z E C
I
- zo)(z - 21).
. .( 5 - zj)l 5 r)
and analytic in the sequence of points { s ~ } ~and ~ Nlet the poles w i of B(z0,. . . , z j , r ) be numbered such that for wi
we have
=
I(Wi
- 20)
0 < w1 Iwp
1 . .
(Wi
- 2j)l
s . .. < r
f in
IV.2.5. The generalized qd-algorithm
237
where the poles are counted with their multiplicities. If the Newton series is ultimately k-normal for some integer k > 0, then the generalized qd-scheme associated with this series has the following properties : a) for each n such that 0 < n 5 k and such that w,-~ < w, < w,+1 where wo = 0 and, if f has exactly k poles, wk+l = r , we have
b) for each ra such that 0 < n 5 k and such that w , < w,+1, we have
For the proof we refer to [9].
238
IV.3.1. The we of Pad4 approzirnants
53. Initial value problems. Consider the following first order ordinary differential equation: dv
- =j(z,y)
dz
for
2
E [a,bj
(4.10.)
with y(a) = yo. When we solve (4.10.) numerically, we do not look for an explicit formula giving y ( z ) as a function of z but we content ourselves with the knowledge of y(zi) a t several points zi in [a, 61. If we subdivide the interval [a,b ] , k
U Izi-1 ,z i ]
[a, 61 =
i= 1
where zi =a
+ ik
i = 0,. . ., k
with
b-a h=--
k
for k
>0
then we can calculate approximations yi+l for y ( z i + ~ )by constructing local approximations for the solution y ( z ) of (4.10.) a t the point 2 ; . We restrict ourselves now t o methods based on the use of nonlinear approximations.
3.1. The use of Pad6 approximantcr. Let us try the following technique. If 8 i ( z ) is the Pad6 approximant of a certain order for y(z) a t z i then we can put Yi+l
=ai(Zi+l)
which is an approximation for y ( z i + l ) . For the calculation of ai(z) we would need the Taylor series expansion of y(z) at z;, in other words
Since the exact value of y(z;) is not known itself, but only approximately by y i , this Taylor series development is not known and hence this technique cannot be applied. However, we can proceed as follows. Consider the power series
IV.3.1. The u8e of Pad6 approzirnants
239
Let r i ( z ) be the Pad6 approximant of order ( m , n )for this power series. If we put z = zi+l, in other words z - x i = h , we obtain 8' = 0,..., k - 1
yi+r = r i ( z i + I )
Hence we can write ui+l
= Y i + hg(zir Yip h)
i = 0,. . .,k - 1
(4.11.)
where g is determined by r i . Such a technique uses only the value of zi and y i t o determine yi+1. Consequently such methods are called one-step methods. Moreover (4.11.) is an explicit method for t he calculation of v i + l . It is called a method of order p if the Taylor series expansion for g(z,y,h) satisfies Y ( z i + l ) - u ( z i ) - hg(ziJ Y ( z i ) ,h) = 0(hp+')
Clearly (4.11.) is a method of order (m + n) if r i ( z ) is a normal Pad6 approximant. The convergence of (4.11.) follows if g(z, y , h) satisfies the conditions of the following classical theorem [30 p. 711.
Theorem 4.8. Let the function g(zyy , h) be continuous and let there exist a constant L such that
g('J YJ
= f(’,
Y)
is a necessary and sufficient condition for the convergence of the method (4.11.), meaning that for fixed z E [ a ,b ] .
2 40
IV.3.1. The we of Podi approzarnants
From the fact that t i is a Pad4 approximant it follows th a t the relation g(z,y, 0) = f(z,y) is always satisfied (see problem (8)). The case n = 0 results in the classical Taylor series method for the problem (4.10.), If we take rn = n = 1 we get (4.12.)
If y(z) is a rational function itself, then using (4.11.) we get the exact solution Yi+l = Y(z;+I)
at least theoretically, if the degrees of numerator and denominator of r i ( z ) are chosen in a n appropriate way. Techniques based on the use of Pad4 approximants can be interesting if we has a large negative real part consider stiff differential equations, i.e. if [ZO]. An example of such a problem is the equation dY = xy dz
(4.13.)
-
with Re(1) large and negative. Since the exact solution of (4.13.) is
we have lim y ( z ) = lim
z-00
2-00
2% =o
and we want our approximations y; to behave in the same way. Dahlquist [15] defined a method to be A-stable if it yields a numerical solution of (4.13.) with Re(1) < 0 which tends to zero as i -+ 00 for any fixed positive h. He also proved t hat there are no A-stable explicit linear one-step methods. Take for instance the method of Euler (rn = 1, n = 0): Y i + l = ~i
+ hf ( z i j y i )
Y i + l = (1
+ hX)Yi
We get
= (1
+ hX)'+1 yo
IV.3.2. The
t~t: of
rational anterpolanta
241
Clearly lim yi+l = yo Iim (I i-
i-+W
00
+ /A)'+' = o
only if
)I
+ hX( < 1
So for large negative X the steplength h has to be intolerably small before acceptable accuracy is obtained. lo practice h is so small t h a t round-off errors and computation time become critical. The problem is t o develop methods t h a t do not restrict the stepsize for slability reasons. If (4.11.) results from the use of t h e Pad6 approximant of order (rn,m), (m, rn 1) or (m,m + 2) then one gets an A-stable method [16]. This can be seen as follows. If f ( z , y) = Xy then r,(z) is the Pad6 approximant for the power series
+
Hence Y i + l = rm,rt(hA) ~i
with h = ~ i + - lz;, where rm+(x) is the Pad6 approximant of order (rn,n)for e 2 . A-stability now follows from the following theorem.
Theorem 4.9. If m = n or m = n - 1 or m = n - 2 then the Pad6 approximant of order (m,n) for ez satisfies
For the proof we refer to [3, 181. 3.2. The uae of rational interpolants.
It is clear t h a t if the interpolation conditions are spread over several points, then the computation of yi+l will need several xi-^ and yi-t(t = 0, 1, . . .). Such methods are called multistep methods. Let r i ( z ) be the rational Hermite iuterpolant of order (rn,n) satisfying
242
IV.8.2. The me of rational interpolants
where (j
+ I)(R + 1) = rn + n + 1
Here
Then an approximation for y(zi+l) can be computed by putting
This is a nonlinear explicit multistep method. A two-step formula is obtained for instance by putting j = 1, m = 2, n = 1, and 8 = 1 and using theorem 3.17.:
where f, = j(x;,y;) and fi-1 = f(z;-1, yi-1). We can also derive implicit methods which require an approximation the caIculation of y;+1 itself, by demanding
v;+l
for
where (j
+ 2)(8 + 1) - 1 = m + n +
For m = 1 = n, j = 0 and
8
1
= 1 we get the formula
(4.15 .)
IV.3.3. Predictor-corrector method8
243
For more information concerning such techniques we refer t o [36]. Remark th a t multistep methods are never selfstarting. Both explicit and implicit ( j + 1)- step methods are of the form Yi+l =
j
C
acyi-c
+ h g ( z i + l , . . ., z i - j l
yi+lJ..
-J
~ i - j lh)
L=O
and they have order p if the Taylor series expansion of g satisfies i
y(zi+t) -
C
aty(zi-t)
- hg(zi+l
J
.
‘1
zi-j] y ( z i + l ) j , ..,y(zi-j), h ) = 0(hp+l)
L= 0
Hence (4.14.) is a third-order method if the starting values are third-order and (4.15.) is second-order. When applied t o a stiff differential equation one should keep in mind that linear multistep methods are not A-stable if their order is greater than two. The following result is helpful if af&’yl is real and negative. We know th a t we can write for problem (4.13.) y(z) = v(2i-j)
e(~-~i-j)X
with Re(X) large and negative. €leiice it is interesting to take a closer look a t rational Hermitc interpolants for exp(z) in some real and negative interpolation points and also in 0. Theorem 4.10.
5 m and 8qi 5 n is such th a t ri( t ) (0) = 1 = exp(*)(O)for 0 5 t 5 rn + n - t2 with e 5 m
If ri(z) = pi(z)/qi(z) with api a)
b) ri(tt) = exp(&) for & then lri(z)l
< 1 if
m
< 0, 1 5 k’ 5 j
5 n and
with &
# &, 1 5 !# k 5 j
z is real and negative.
For the proof we refer to [32]. 3.3. Predict0 r-correct or methode. For a solution of the initial value problem (4.10.) we have y’(z) = f ( z , y(z)) for every z in [a, b]
244
IV.S.4. Numeracal results
If we integrate this equation on the interval [zi-i, z i + t ] with j , 8 2 0 we get PZitl
Now f can be replaced by an interpolating function, through the points j(zi-1 ,yj-l)), . . ., which is easily integrated. ( z 8 ,f(z;, yi)), If t! = 1 we get p r e d i c t o r - m e t h o d s because they are explicit. If 8 = 0 we get corrector-methods because the value of y; is needed for the computation of y;. These implicit formulas can be used t o update an estimate of yfz;)iteratively. W h e n f is replaced by an interpolating polynomial we get the well-known methods of Adarns-Bashforth ( j = 0 and !. = 1) and Adams-Moulton ( 3 = 1 and 6 = 0). When f is replaced by an intc3rpolating rational function we can get nonliriear formulas of predictor or corrector type. Sincc the rational interpolant must be integrated, it is not recommendable to choose rational functions with a denominator of high degree.
3.4. Nu me rical r es ult 8 . Let u s compare two Taylor series methods ( m= 2 and m = 3) with the explicit method (4.12.) for t h e solution of the equation y‘ = I
+ y2
for z
2o
Y(0) = 1
The theoretical solution is g = tg(z -t 7r/4). We take the steplength h = 0.05. As can be seen in table 4.19. the second order rational method gives even better results t h a n t h e third order Taylor series method, a fact which can be explained by the singularity of the solution y(z) at z = 7 . To illustrate A-stability we will compare Euler’s method ( m= 1,n = 0) with the formulas (4.12.))(4.14.) and (4.15.) for the equation y’ = -25y
for z E [0,I]
Y(0) = 1
The solution is known to be y = exp(-25z). For the results in table 4.20. we chose the steplength h = 0.1. As expected Eiiler’s solution blows up while formulas based on the use of a “diagonal” entry of the Pad6 table or the rational Hermite interpolation table for t h e exponential decay quite rapidly. Similar results would have been obtained if “superdiagonal” entries of the Pad6 table were used. Remark that formula
Table 4.19. exact solution
+7)
Taylor series
Pad6 approximant
Taylor series
m=3
i
zi
1
0.05
0.1 10536D+01
0.1 105OOD+Ol
0.110526D+01
0.110533D+01
2
0.10
0.122305D+01
O.l22210D+Ol
0.122284D-kOl
0.122299D)+Ol
3
0.15
0.135600D-kOl
0.135449D+Ol
0.135573D+01
0.135598D)+01
4
0.20
0.150850D+01
0.1 5 0 5 8 2 D t 0 1
0.150795D+01
0.150831D+Ol
5
0.25
0.1 6858oD+oi
0.168150D-kOl
O.l68500D+Ol
0.188547D+01
6
0.30
0.188577D$-01
O.l88886D+Ol
0.188462D-kOl
0.189522D-t-01
7
0.35
0.214875D-kOl
0.213884D-kOl
0.214811D+Ol
0.214883D-kOl
8
0.40
0.246496D-kOl
0.244751D-kOl
0.24626lD+Ol
0.246335D-kOl
9
0.45
0.286888D+01
0.28398OD+Ol
0.286543D-I-01
0.286584D-I-01
10
0.50
0.340822D+Ol
0.335737D-t-01
0.340298D-t-01
0.340248D-kO 1
11
0.55
0.4 16936D4-01
0.407388D-kOl
0.416097D+Ol
0.415703D+01
12
0.60
0.533186D-t-01
0.513307D+Ol
0.531720D+Ol
0.5301 31D-t-01
13
0.65
0.7 3 4 0 4 4 D t 0 1
0.685144D-t-01
0.731087D+01
0.724568D-kOl
14
0.70
0.1 16814D+02
0.1 00697D-kO2
0.1160 1QDi-02
0.112431D-kO2
15
0.75
0.282383D-kO2
0 .I 77676D-kO2
0.277486D-i-02
0.232131D-kO2
tg(zi
m=2
m=l=n
*CD e
5 e
Table 4.20.
exact solution exp(-25si)
i 1
0.1
0.820850D-01
2
0.2
0.6737851)-02
3
0.3
0.553084D-03
4
0.4
0.453QQQD-04
5
0.5
0.372665D-05
6
0.6
0.305902D-06
7
0.7
0.251 100D-07
8
0.8
0.2061 15D-08
Q
0.Q
0.16QlQOD-OB
10
1.0
0.13887QD--10
m=l=n explicit one-step
Euler -0.150OOOD+01 0.225000D+01 -0.337500D+01 0.506250D-l-01 -0.758375D+01 0.113Q06D+02 -0.17085QD)+02 0.25628OD-l-02 -0.384434D)+02 0.576650D4-02
1
-0.111111D+OO 0.123457D-01 -0.137174D-02 0.152416D-03 -0.189351 D-04 0.188168D-05 -0.200075D-06 0.2323061)-07 -0.258117D-08 0.286797D-00
rn = 2,n = I
explicit multisteD
rn=l=n
implicit one-step
0.820850D-01
0.123047D+00
0.840728D-01
0.151 407D-01
-0.542641D-01
0.186302D-02
-0.4Q474OD-kOO
0.22023QD-03
-0.252173D+00
0.282073D-04
0.314055D-kOO
0.347083D-05
0.102175D+O1
0.427077D-06
0.16Ql37D-kOO
0.525507D-07
-0.lQl5Q6D-I-00
0.646622D-08
-0.6Q8115D-kOO
0.7Q5651D-OQ
IV.3.5. Systems of first order ordinary differential equatiow
247
(4.14.) based on a "subdiagonal" rational approximation is surely not producing an A-stable method. To obtain the results of table 4.20. by means of (4.14.) a second starting value yl was necessary. We took yl = exp(-2.5) = y(z1). For (4.15.) the expression f(zi,y;) = -25y; was substituted and y;+l was solved
from the quadratic equation. All these schemes can be coupled to mesh refinement and the use of extrapolation methods. If an asymptotic error expansion of y(zi) in powers of h exists, then the convergence of the sequence of approximations for y(z;), witb z i fixed, obtained by letting the stepsize decrease, ran be accelerated by the use of techniques described in section 1 [23]. 3.5. Systems of first order ordinary differential equations.
Nonlinear techniques can also be used to solve a system of first order ordinary differential equations
dzj = f j ( X , 2 1 , 2 2 , . . . , Z t ) dX
where the values z,,! and the functions approaches are possible. If we introduce vectors
fj
j = 1, ..., k are given for j = I , . . .,k. Several
then one method is t o approximate the solution componentwise using similar techniques as in the preceding sections. So for instance (4.12.) becomes
Y ( z ; + ~N) Y;:+l where
and
= Y;
Y;) ] + h [*F(zi,2F2(Zi, Y;)- hF'(xil yi)
(4.16.)
248
IV.3.5. Systems of first order ordinary differential equation8
and the addition and multiplication of vectors is performed componentwise. In other words (4.16.) is equivalent with
For more information on such techniques we refer to [56,35, 381. Another approach is not based on componentwise approximation of the solution vector Y ( s )but is more vectorial in nature. Examples of such methods and a discussion of their properties is given in [57, 7 , 27, 131. The nonlinear techniques introduced here c a n also be used to solve higher order ordinary differential equations and boundary value problems because these can be rewritten as systems of first order ordinary differential equations. Again same of the nonlinear techniques prove to be especially useful if we are dealing with stiff problems. A system of differential equations
dY -= F ( z , Y ) dx Y ( a ) = Yo with
is called stiff if the matrix
IV.3.5. System8 of first order ordinary differential equations
249
has eigenvalues with small and large negative real part.
{
Consider for example
+
= 9 9 8 ( ~x )~ 199822 (z)
Zl(0j = 1
- - - -99921(2) - 19992z(x)
22(0)
=0
The solution is
so that again
lim z I ( x ) = 0 = lim
z 4 m
DO’%
z2(x)
where both z1 and a2 contain fast and slow decaying components. For a discussion of stiff problems we refer to [20 pp.209-2221.
IV.4. Numerical integration
250
§4. Numerical integration.
s,
b
Consider I = j ( z ) d z . Many methods to calculate approximate values for I are based on replacing j by a function which can easily be integrated. The classical Newton-Cotes formulas are obtained in this way: j is replaced by an interpolating polynomial and hence I is approximated by a linear combination of function values. In some cases the values of the derivatives of j ( x ) are also taken into consideration and then linear combinations of the values of f ( z ) and its derivatives a t certain points are formed t o approximate the value I of the integral. This is for instance the case if polynomial Hermite interpolation is used. In many cases the linear methods for approximating I give good results. There are however situations, for example if f has singularities, for which linear methods are unsatisfactory. So one could t r y t o replace f by a rational function r and consider
[
r(x)dz
as an approximation for I . But rational functions are not that easy to integrate unless the poles of r are known and the partial fraction decomposition can be formed. Hence we use another technique. Let us put
Then
1 = y(b) If f is Riemann integrable on [a,b ] , then y is continuous on [a, b ] . If f is continuous on [a, 61, then y is differentiable on [a,b] with y’(z) = j ( z ) and y(a) = 0
So I can be considered as the solution of a n initial value problem and hence the techniques from the previous section can be used. We group them in different categories.
I V . 4 . l . Methoda w i n g Pad6 approximantd
251
4.1. Methods using Pad6 approzirnants. Let us partition the interval [a, 61 with steplength h = ( b - a ) / k and write Z; = a
+ ih
a = 0, ...,k
and (4.17.)
s,”’
where y ; approximates y ( z i ) = f(t)dt. If r; is the Pad6 approximant of order ( m ,n ) for t l , ; ( h )then we can put
and consider y k as an approximation for I . In this way Yi+l
which means
= yi
+ hg(zi, h)
i = 0, ...,k - 1
(4.18.)
l;’’
f ( t ) d t N hg(z;,h)
If m = 1 = ra we can easily read from (4.12.) t h a t (4.18.) results in (4.19.) Formulas like (4.18.) use derivatives of f(z) and are nonlinear if n > 0. From the previous section we know that (4.18.) is exact, in other words that yk = I , if y ( z ) is a rational function with numerator of degree m and denominator of degree ra. For n = 0 formula (4.18.) is exact if y ( z ) is a polynomial of degree m, i.e. f ( z ) is a polynomial of degree m - 1. The obtained integration rule is then said t o be of order m - 1. The convergence of formula (4.18.) is described in the following theorem which is only a reformulation of theorem 4.8.
252
IV.4.2. Method8 wing rational interpolanta
Theorem 4.11. Let y; be defined by (4.18.). Then lim h-0
yi(h) = y(z) for fixed z E [a, b]
====i(h)
if and only if g(z,O) = f(z). Instead of (4.17.) one can also write t l , i ( h )= T/i
+ lata,i(h)
with
and compute Pad6 approximants r , of a certain order for tz,;(h).Let us take
m = 1 = n . If we define
Yi+i
= yi
+ hr;(h)
then we get
4.2. Method8 using rational interpolanta.
If we again proceed as in the section on initial value problems we can construct nonlinear methods using information in more than one point. Since these methods are not self starting but need more than one starting value their use is somewhat limited and rather unpractical. An example of such a procedure is the following. Let ti(.) be the rational Hermite interpolant of order ( 2 , l ) satisfying
and let
IV.4.3. Methods without the evaluation of derivatives
253
Then we know from formula (4.14.) t h a t
with f; = f(zi) and f i - 1 = f ( z i - I ) . Often rational interpolants are preferred t o Pad6 approximants for the solution of numerical problems because the use of derivatives of f is avoided. As mentioned, a drawback here is the necessity of more starting values. Another way t o eliminate the use of derivatives, now without the need of more starting values, is the following.
4.3. Methods without the evaluation of derivatives. One can replace derivatives of j ( z ) in formula (4.18.) by linear combinations of function values of f ( z ) without disturbing the order of the integration rule. To illustrate this procedure we consider the case m = 1 = n. Then 1591
We will compute constants a,p and 7 such t h a t for
we have gfZ,
h) - t ( s ,h) = O(hrn+%) = O(h2)
For ~ ; + i= ~i
+ h t ( z ; ,h)
this would imply y(si+l) - Yi+l = ~(hm+m+l)= O ( P )
Condition (4.20.) is satisfied when a+B=2 B7 = -1
(4.20.)
254
IV.4.4. Numerical resulta for singular integranda
In other words, for 7 # 0, 27 a=---
p = -1
+1 7
-.
7
so
For 7 = 1 we get the integration rule
In this way we approximate (4.21.)
4.4. Numerical results for singular integranda We will now especially be interested in integrands regular in [a,b] but with a singularity close t o the interval of integration and on the other hand in integrands singular in Q or b . The problem of integrating a function with several singularities within [a, 61 can always be reduced t o the computation of a sum of integrals with endpoint singularities. If f ( z ) is singular in 6 then the value of the integral is defined by
and is assumed t o exist. We shall compare formula (4.19.) with Simpson’s rule ( m = 2, n = 0 ) and with a (2k)-point Gaussian quadrature rule that isolates the singularity in the weight function. If f(t) can be written as w ( t ) h ( t ) where w ( t ) contains the endpoint singularity of f(t) and h(t) is regular then the approximation
I
wlh(tl) + . . .
+ 402kh(t2k)
does not involve function evaluations in singular points. We use a (2t)-point formula because on [zi,2;+1] for i = 0 , . . ., k - 1 both (4.19.) and Simpson’s rule
IV.4.4. Numerical reaulta for aingular integranda
255
need two function evaluations. Since / is singular in b = zk, we take f ( z k ) = 0 in Simpson's rule which means that the singularity is ignored. In (4.19.) the singularity of f is no problem since f and f’ are only evaluated in 2 0 , . . ., xk--1 [58]. Our first numerical example is
d t = 3.04964677.. with a singularity in
t = In 3 = I .09861228. . . and the second example is
+
2(1_ _ ~ t)sint _ _ _ _ _cost d t = 2
We shall also compare the different integration rules for the calculation of e'dt = 1.71828182. . .
which has a smooth integrand. Because of the second example the weight function w ( t ) in the Gaussian quadrature rule was taken t o be W ( t )=
1
J1-t
.___
All the computations were performed in double precision arithmetic and the double precision values for the weights w i (i = 1 , . . .,2k) and the nodes ti (d = 1 , . . .,2k) were taken from [I pp. 916-9191. Remark t ha t the nonlinear techniques behave better than the linear techniques in case of singular integrands. However, for smooth integrands such as in table 4.23., the classical linear methods give better results th a n the nonlinear techniques. Also, if the singularity can be isolated in the weight function such as in table 4.22., Gaussian quadrature rules are very accurate. In general, little accuracy is gained by using nonlinear techniques if other methods are available for the type of integrand considered [51].
256
IV.4.4. Numerical reeults f o r singular integranda
As in the previous section all these schemes can be coupled to mesh refinement and extrapolation.
Table 4.21. 1
I=J
et
o (3 - et)2
k = 16
dt = 3.04964677.. .
Gaussian
Simpson's
formula
quadrature
rule
(4.19 .)
3.1734660 3.0773081 3.0564777
3.28067 65 3.0841993 3.0531052
3.1006213 3.0573798 3.0510553
Table 4.22.
I = J 2(1 - t ) sin t -tc o s t dt = 2 0 G - j
k=4 k=8 k = 16
Gaussian
Simpson's
formula
quadrature
rule
(4.19.)
2.0000000 2.0000000 2.0000000
1.7504074 1 A267581 1.8786406
1.9033557 1.8797832 1.go48831
Table 4.23. 1
I = 1etdt = 1.71828182.. . 0
k=4 k=8 k = 18
Gaussian
Simpson's
formula
quadrature
rule
(4.18.)
1.7205746 1.7204038 1.7188196
1.7182842 1.7182820 1.7182818
1.7284991 1.7206677 1.7188592
IV.5. Partial diflerential equatiom
257
§5. Partial differential equations. Nonlinear techniques are not frequently used for the solution of partial differential equations. We will describe here a method based on the use of Pad6 approximants t o solve the heat conduction equation which is a linear problem. For other illustrations we refer to [54, 48, 17, 22, 281. All these techniques first discretize the problem so that thr: original partial differential equation is replaced by a system of equations which is nonlinear if the partial differential equation is. Another type of techniques which we do not consider here are methods which do not discretize the original problem but solve it iteratively by means of a procedure in which subsequent iteration steps are differentiable functions [13]. Linear techniques of this type arc recommended for linear problems and nonlinear techniques can be used for nonlinear partial differential equations. Let us now concentrate on the heat conduction equation. Suppose we want t o find a solution u ( z , t )of the linear problem
a+, t ) - a%+, t ) ~
at
.~
82
~
a
t>O
(4.22 .)
with u(2,O) = u ( 2 )
u(0,t ) = a
4 4 t)= P The domain [a, 61 with
x
[0,oo) is replaced by a rectangular grid of points (z;, ti)
zi = a + ih ti = j A t
i = 0,.. .,k + 1
j = 0, 1 , 2 , . . .
where
and At is the discretization step for the time variable. We first deal with the discretization in the space variable z and introduce u;(t) = u ( z ; , t ) for
t20
Then using central differences, (4.22.) can for instance be reduced t o
258
I V.5. Part id digere nt ial
equ at ion8
for t > 0 and i = 1, . . . , k. In general we c a n write
(4.23.)
where A is a real symmetric positive definite k x k matrix and depends on the chosen approximation for the operator d 2 / d s 2 . If we introduce the notations
then the exact solution of (4.23.) is U ( t ) = eWtAV where e-tA is defined by
Using the discretization in the time variable t we can also write (4.24.)
I V.5. Partial d ifler e nt ial equ at i o ns
259
If rm,.,(t)= p ( t ) / q ( t ) is the Pad6 approximant of order (n,n) for e-t then (4.24.) can be approximated by
U ( t + A t ) = [q(At.A)]-' [P(At.A)] U ( t )
(4.25 .)
Varga proved t ha t for n 2 m this is an unconditionally stable method [ 5 5 ] , meaning t ha t initial rounding errors remain within reasonable bounds as the computation proceeds independent of the stepsize used. If rn = 1 and n = 0 then (4.25.) means
U(t + At) = ( I
- AtA)U(t)
or equivalently if A = (&.!)kXk
ui(ti+l) = u i ( t j ) - ~t
k
C aitut(tj)
.!=I
which is the well-known explicit method to solve (4.22.). The solution a t level t = ti+l is determined from the solution at t = t , , For nt = 1 = n we obtain from (4.25.)
or equivalently
which is the method of Crank-Nicholson [50]. The operator a2/ax2is replaced by the mean of a n approximation for the partial second derivative a t level t = and the same approximation at t = t i .
IV.6.1. Kernel8 of finite rank
260 56. Integral equations.
As in the previous section we shall discuss linear equations for which the use of nonlinear methods is recommendable. Those interested in nonlinear integral equations are referred to [lo, 121 where methods are indicated for their solution. If the integral equation is rewritten as a differential equation then techniques developed for the solution of initial value problems can also be used. We restrict ourselves t o the discussion of an inhomogeneous Fredholm integral equation of the second kind (the unknown function f appears once outside the integral sign and once behind it):
f ( 4-
J'." K ( z ,y)f(v)dy
= g(z)
xE
1% bl
(4.28.)
Here the kernel K ( z , y ) and the inhomogeneous right hand side g(z) are given real-valued continuous functions. Fredholm equations reduce to Volterra integral equations if the kernel K ( z , y ) vanishes for g > z which produces a variable integration limit.
6.1.K e r n e l 8 of finite rank. Formally the solution of (4.26.) can be written as a series. P u t
and
If we define
and
then (4.26.) reduces to
IV.6.1. Kernel8 of finite rank
261
The series (4.27.) which is a power series in X, is called the Neumann series of the equation (4.26.). Convergence of the Neurnann series for certain values of X depends on the properties of the kernel K ( z ,y). If K ( z ,y ) is bounded by IWz,y)l
for (z,Y) E
b,b]
X [a, b ]
then clearly the series (4.27.) converges uniformly t o f(s)in [a, b] if 1
1x1 < ---___ M ( b - a) If Pad6 approximants in the variable X are constructed for (4.27.) then they may have a larger convergence region than the series itself. Especially interesting is the case t ha t the kernel is degenerate, in other words k
~ ( zY), =
C Xi(Z)Yi(y)
i= 1
with ( X i } and {Y;} each linearly independent sets of functions. Such a kernel is also said to be of finite rank k. Let us try to determine f(z) in this case. If we put
then (4.26.) can be written as (4.28.) Multiply (4.28.) by Yj and integrate to get
or equivalently
IV.6.2. Completely continuoua kernela
262
(4.29 .)
which is the determinant of the coefficient matrix of system (4.29.), then D(X) is a polynomial in X of degree at most k. In case D(X) # 0 the solution of (4.29.) is given by
with
Dij
the minor of the ( i , j ) t h element in D(X). SO (4.28.) can be written as
which is a rational function in X of degree a t most k both in numerator and denominator. Since the series development of (4.30.) coincides with the Neumann series (4.27.) we know that the solution f(z) is equal to the Pad6 approximant of order (k,k) for the Neumann series. So in this case the Pad6 approximant is the exact sum of the series because the sum is a rational function [2 p. 1761. 6.8. Completely continuoua kernels.
The equation (4.26.) can be rewritten
EUI
(I-XK)f=g where the linear operators I and K are defined by
(4.31.)
263
IV.6.2. Completely continuoua kernela
We shall consider square-integrable functions f with
Suppose now t ha t { f;}iEN is a bounded converging sequence of functions. Then the operator K is said t o be completely continuous if for all bounded { f ; ) , ~ the sequence { K fi}iENcontains s subsequence converging t o some function h ( s ) with
- h(2)l2ds -+ 0 as i -+
llKfi - hll = \i[lKfi(Z)
00
A basic property of completely continuous transformations is th a t they can be uniformly approximated by transformations of finite rank. Thus there is a n infinite sequence { K ; } ; Lof~ kernels of finite rank i such th a t
where lim
i-+w
t;
=0
When a completely continuous kernel K is replaced by Ki then the solution f; of
( I - Mi)!; = g is given by the Pad4 approximant ri,i of order (i,;), as explained in the previous section. In case K is completely contiuuous the exact solution f of (4.31.)is a meromorphic function of X [31 p. 311 arid
lim f, = f
i-
00
2 64
Iv.6.2. Completely continuous kernels
in any compact set of the X-plane except at the finite number of X-poles of f [8].Since
we consequently have
lim
ri,i
1-m
=f
in any compact set of the X-plane except at the finite number of X-poles of f.Thus the solution f is the limit of a sequence of diagonal Pad6 approximants. This result is not very useful since each approximant in the sequence is derived from a different Neumann series, with kernel K ; and not with kernel K . However a similar result exists for the sequence r i , i derived from the Neumann series (4.32.) i= 1
where the operator K' is defined by
Theorem 4.12. If ri,i is the Pad4 approximant of order (d, i) t o the Neumann series (4.32.) of the integral equation (4.28.) with completely continuous kernel K then the solution f is given by lim
i-oc
ri,i
=f
in any compact set of the X-plane except a t the finite number of poles of f and at limit points of poles of ri,i as i -+ 00. More information on this subject can be found in [2 pp. 178-1821.
I V . Problems
265
Problems. Show t ha t with
EF’= an for n 2 0,
which coincides with Aitken’s A2-process to accelerate the convergence of a sequence. Give a n algorithm similar to the one in section 1.4. for the calculation of t,,i, but now based on the use of inverse differences instead of reciprocal differences, Show that if the algorithm of Bulirsch and Stoer is used with the interpolation points zi (lim,+m s; = 0) and the p-algorithm is used with the interpolation points z: = 1,s; (1imi4- 2;. = oo),then t i ; = B~( 0 ) . Compare the amount of additions and multiplications performed in the c-algorithm and the qd-algorithm when used for convergence acceleration. Derive the formulas (4.3.) and (4.4.) using inverse interpolation instead of direct interpolation. Derive the formulas (4.8.) and (4.9.) based on the use of rational interpolants. in section 2.3. for successive values of Organize the computation of d,+l d such t ha t a mininum number of operations is involved. Prove that one-step explicit methods for the solution of initial value problems based on the use of Pad6 approximants are convergent if g(z, g, h) given by (4.11.) is continuous and satisfies
IV. Problem8
266
(9)
Check the formulas (4.12.) and (4.14.).
(10)
Write down formula (4.16.) for the solution of
- 99821 (2)+ 199822( Z) - - - -9992+) (11)
- 199922(.)
ZI(0) = 1 ZZ(0) = 0
Construct a nonlinear numerical integration rule based on the use of Pad4 approximants of order ( 2 , l ) . Afterwards eliminate the use of derivatives as explained in section 4.3.
I V . Remnrka
267
Remarks. Nonlinear methods can also be used for the solution of other numerical problems. We refer for instance to [40] where the solution of linear systems of equations is treated, t o [19] for analytic continuation, to [34] for numerical differentiation and to (421 for the inversion of Laplace transforms.
An important link with the theory of numerical linear algebra is through QR-factorization. Rutishauser [47] proved th a t the determination of the eigenvalues of a square malrix A can be reduced to the determination of the poles of a rational functiom .f built from th e given matrix. In this way decomposition techniques for A to compute its eigenvalues are related t o the qd-algorithm when used to compute poles of meromorphic functions. Univariate and multivariate continued fractions and rational functions are also often used t o approximate given functions. For univariate examples we refer t o (33, 421. The bivariate Beta function is a popular multivariate example because it has numerous singularities in a quite regular pattern. For numerical results we refer to [89, 26, 121.
As a conclusion we can say t h a t every linear method has its nonlinear analogue. In case linear methods are inaccurate or divergent, it is recommendable t o use a similar nonlinear technique. The price we have to pay for the ability of the nonlinear method to cope with the singularities is the programming difficulty to avoid division by small numbers.
2 68 References.
1
1) Abramowite M. and Steguo Handbook o Mathematical functions. Dover publications, New York, 1968.
[
21 Baker G. and Gamrnel J. The Pad6 approximant in theoretical physics. Academic Press, New York, 1970.
[
31 33aker G. and Graves-Morris P. Pad4 approximants: basic theory. Encyclopedia of Mathematics and its applications: vol 13. AddisonWesley, Reading, 1981.
\
41 Brezinski C. Algorithnies d’accklkration de la convergence. Editions Tecbnip, Paris, 1978.
[
51 Rrezinski C. Application de 1’6-algorithme B la r6solution des syst6mes non linhaires. C. R. Acad. Sci. Paris Sbr. A 271, 1970, 1174-1177.
[
61 Bulirsch R . and Stoer J. Fehlerabschatzungen und Extrapolation mit rationalen Funktionen bei Verfahren vom Richardson-Typus. Numer. Math. 6, 1964, 413-427.
[
On the stability of rational Runge 71 Calvo M. and Mar Quemada M. Kutta methods. J. Comput. Appl. Math. 8, 1982, 289-292.
[
81 Chisholm J. Solution of linear integral equations using Pad6 Approximants. J. Math. Phys. 4, 1963, 1506-1510.
[
91 Claessens G. Convergence of multipoint Pad6 approximants. Report 77-26, University of Antwerp, 1977.
[ lo] Clarysse 7 .
Rational predictor-corrector methods for nonlinear Volterra integral equations of the second kind. In [60], 278-294.
[ 111 Cuyt A.
Accelerating the convergence of a table with multiple entry. Numer. Math. 41, 1983, 281-286.
[ 121 Cuyt A.
Pad6 approximants for operators: theory and applications. Lecture Notes in Mathematics 1065, Springer Verlag, Berlin, 1984.
[ 131 Cuyt A.
Pad6 approximants in operator theory for the solution of nonlinear differential and integral equations. Comput. Math. Appl. 6, 1982, 445-466.
Iv. Reference8
269
Abstract Pad6 approximants for 141 Cuyt A. and Van der Cruyssen 1’. the solution of a system of nonlinear equations. Comput. Math. Appl. 9, 1983, 617-624.
[ 151 Dahiquist G.
A special stability problem for linear multistep methods. BIT 3, 1963, 27-43.
[ 161 E61e B.
A-stable methods and Pad4 approximations t o the exponential. SIAM J. Math. Anal. 4, 1973, 671-680.
[ 171 Fairweather G.
A note on the eficient implementation of certain Pad6 methods for linear parabolic problems. BIT 18, 1978, 100-108.
[ 181 Frame J .
The solution of equations by continued fractions. Amer. Math. Monthly 60, 1953, 293-305.
[ 191 Gammei J.
Continuation of functions beyond natural boundaries. Rocky Mountain J. Math. 4, 1974, 203-206.
[ 201 Gear C .
Numerical initial value problems in ordinary differential equations. Prentice-Hall Inc., New Yersey, 1971.
[ 211 Genz A .
The approximate calculation of multidimensional integrals using extrapolation methods. Ph. 1). in Appl. Math., University of Kent, 1975.
[ 221 Gerber P. and Miranker W .
Nonlinear difference schemes for linear partial differential equations. Computing 11, 1973, 197-212.
[ 231 Gragg W.
On extrapolation algorithms for ordinary initial value problems. SLAM J. Numer. Anal. 2, 1965, 384-403.
[ 241 Gragg W. and Stewart G .
A stable variant of the secant method for solving nonlinear equations. SIAM J. Numer. Anal. 13, 1976, 889-903.
[ 251 Graves-Morris P .
Pad6 approximants and their applications. Academic
Press, London, 1973.
[ 261 Graves-Morris P. , Hughes Jones R . and Makiuson G . The calculation of some rational approximants in two variables. J. Inst. Math. Appl. 13, 1974, 311-320.
[ 271 Hairer E.
Nonlinear stability of RAT, an explicit rational Runge-Kutta method. BIT 19, 1979, 540-542.
IV. References
270
1
Pad6 approximants, fractional step methods 281 Hall C. and Porsching T. and Navier-Stokes discretizations. SLAM J . Numer. Anal. 17, 1980, 840851.
[ 291 Henrici P.
Applied and computational complex analysis: vol. 1. John Wiley, New York, 1974.
[ 301 Henrici P.
Discrete variable methods in ordinary differential equations. .John Wiley, New York, 1962.
[ 311 Hoheisef G.
Integral equations. Ungar, New York, 1968.
[ 32) Iserles A.
On the generalized Pad6 approximations to the exponential function. SLAM J. Math. Anal. 16, 1979, 631-636.
[ 33J Kogbethnts E .
Generation of elementary functions. In [46], 7-35.
[ 341 Kopal Z.
Operational methods in numerical analysis based on rational approximation. In 1371, 25-43.
[ 351 Lam bert J.
Computational methods in ordinary differential equations. John Wiley, London, 1973.
A method for the numerical solution of [ 361 Larnbert J. and S6aw B. g' = f ( z , y) based on a self-adjusting non-polynomial interpolant. Math. Comp. 20, 1966, 11-20.
[ 371 Langer R.
On numerical approximation. University of Wisconsin Press, Madison, 1959.
[ 381 Lee D.and Preiser S.
A class of nonlinear multistep A-stable numerical methods for solving stiff differential equations. Internat. J. Comput. Math. 4, 1978, 43-52.
I 391 Levin D.
On accelerating the convergence of infinite double series and integrals. Math. Comp. 35, 1980, 1331-1345.
[ 401 Lindskog G.
The continued fraction methods for the solution of systems of linear equations. BIT 22, 1982, 519-527.
[ 411 Longman I.
Use of Pad4 table for approximate Laplace transform inversion. In [25], 131-134.
The special functions and their approximations. Academic 421 Luke Y . Press, New York, 1969.
IV. References
271
[ 431 Merz G .
PadCsche Naherungsbruche und Iterationsverfahren hoherer Ordnung. Computing 3, 1968, 165-183.
[ 441 Nourein M.
Root determination by use of Pad6 approximants. BIT 16, 1976, 291- 297.
[ 451 Ortega J. and RheinboIdt W.
Iterative solution of nonlinear equations in several variables. Academic Press, New York, 1970.
[ 461 Ralston A. and Wilf S.
Mathematical methods for digital computers. John Wiley, New York, 1960.
Der Quotienten-Differenzen AIgorithmus. [ 471 Rutishauser H. Math. Phys. 5, 1954, 233-251.
z. h g e w .
On time-discretiEations for linear time-dependent partial differential equations. University of Manchester, Numer. Anal. Report 5, 1974.
[ 481 Siemieniuch J . and Gladwell I .
Numerical comparison of nonlinear convergence accelerators. Math. Comp. 38, 1982, 481-499.
[ 491 Smith D. and Ford W. [ 501 Smith G.
Numerical solution of partial differential equations. Oxford University Press, London, 1975.
[ 511 Squire W.
In defense of linear quadrature rules. Comput. Math. Appl. 7, 1981, 147-149.
[ 521 Tornheim L.
Convergence of multipoint iterative metho&. J . Assoc. Comput. Mach. 11, 1964, 210-220.
[ 531 n a u b J .
Iterative methods for the solution of equations. Prentice-Hall Inc., New York, 1964.
[ 541 Varga R.
Matrix iterative analysis. Prentice-Hall Inc., Englewood-
Cliffs, 1962.
[ 551 Varga R .
On high order stable implicit methods for solving parabolic partial differential equations. J . Math. Phys. 40, 1961, 220-231.
[ 561 Wambecq A.
Nonlinear methods in solving ordinary differential equations. J. Comput. Appl. Math. l , 1976, 27-33.
Iv. References
272
[ 571 Wambecq A.
Rational Runge-Kutta methods for solving systems of ordinary differential equations. Computing 20, 1978, 333-342.
[ 581 Werner H . and Wuytack L.
Nonlinear quadrature rules in the presence of a singularity. Comput. Math. Appl. 4, 1978, 237-245.
[ 591 Wuytack L.
Numerical integration by using nonlinear techniques. J. Comput. Appl. Math. 1, 1975, 267-272.
[ 601 Wuytack L.
Pad6 approximation and its applications. Lecture Notes in Mathematics 765, Springer, Berlin, 1979.
1
611 Wynn P . Acceleration techniques for iterated vector and matrix problems. Math. Comp. 16, 1962, 301-322.
[ 621 Wynn P. 175-195.
Singular rules for certain non-linear algorithms.
BIT 3, 1963,