Post on 06-Nov-2015
description
transcript
Numerische Mathematik 6, 413 -- 427 (1964)
Feh lerabsch~i tzungen und Ext rapo la t ion
mi t ra t iona len Funkt ionen
be i Ver fahren vom R ichardson-Typus
Yon
ROLAND BULIRSCH und JOSEF STOER
1.
Sei T(h) die zu einem Diskretisierungsparameter h gewonnene numerische N/iherung eines exakten Problems, definiert durch lim T(h)= T(0). Von RI-
h--~0 CHARDSON [8] stammt folgende Idee einer Verbesserung der T(h): Man berechne T(h 3 far verschiedene hi, i=0 . . . . . m, lege durch die Stiitzpunkte (hi, T(hi) ) ein interpolierendes Polynom Tin(h) und nehme T~,(0) als N/iherungswert fiir T(0). LAURENT [7] hat ktirzlich in einem allgemeinen Rahmen gezeigt, dab unter geeigneten Voraussetzungen Tin(0) mit wachsendem m gegen T(0) kon- vergiert.
RICHARDSONs Methode, verbunden mit einem von NEVILLE bzw. AITKEN stammenden Interpolationsalgorithmus zur Ermittlung von Tm (0), liefert im all- gemeinen mit geringem Aufwand sehr gute numerische Resultate. Beispiel daf~lr ist das bekannte, nach ROMBERG [91 benannte Quadraturverfahren; vgl. dazu die Arbeit [1] von BAUER, RUTISHAUSER und STIEFEL und die Arbeit [3]. Weitere Anwendungen finden sich in Arbeiten yon RUTISHAUSER [103, BOLTON and SCOINS [23 u.a. Siehe dazu auch die Arbeiten [13, 14, 15, 1611.
Voraussetzung ftir die numerische Brauchbarkeit der Extrapolationsmethode ist atlerdings die Existenz einer Entwicklung der Form
(l) T(h) = ~0 + ~1 h~' + ... + ~k h~ + R~+I (h) h~
mit [R~ fiir alle h>0, % . . . . . z~ unabh/tngig von h und 0
4t4 ROLAND BULIRSCH und JOSEF STOER:
Formel angegeben werden. Beispiele zeigen die ~3berlegenheit dieses Verfahrens, das, obwohl nur geringftigig komplizierter, in allen untersuchten Fiillen nicht schlechter, meistens sogar erheblich besser konvergierte als entsprechende Poly- nom-Verfahren.
2.
Zun~chst sollen lineare Extrapolationsoperatoren A~ mit den Eigenschaften
i+m a) A~ T = T(m i) -- ~ c(O T(hj)
9 . ~ml
(2) b) A~ t = 1 '='
c) A~h~'=O, i 1 . . . . . m
rekursiv konstruiert werden. Dabei gelte ho>ha:>-. .>0. Definiert man den Ausdruck (in Spezialffillen das Polynom)
i ' h J
durch die Forderung
T,,,(hi)= T(hi), j= i , i + i . . . . . i +m,
so ist die Bestimmung yon A~T gleichwertig der Berechnung von T~(0). Ist n~imlich T,,(h) ein solcher Ausdruck, so gilt wegen (2)
AmT, ,=A, ,T=ao=T, , (O) .
Zur Konstruktion der A~m lRl3t sich also die Theorie der Interpolation heranziehen. Fiir den Fall Yi=/'Y liefern die Nevilleschen Interpolationsformeln, angewandt auf ein Polynom in h v (s. [5], ~6]), die L6sung
(3)
A~ T---- T0(0= T(hi),
A~ T---- T~ ~) h~T'(i+i)-- h~ T'(i) = ~ m- -1 i+m m--1 h~- h L r~
__ T ( i + I )
(2,o)'-' Ordnet man die T~ ) in dem Schema an
T(%) = To (~ T:( ~ )
T(h:) = 1"(o :1 T(2 ~ T1 (1)
T(h2) = To(~) T~(:) TI(2)
T(ha) = T~ ~)
(4)
, m~ t .
7~ ~ )
so lassen sich die T~ ) ausgehend vonder ersten Spalte rekursiv berechnen. Fiir die c(~)/ aus G1. (2) liefert die Interpolationsformel yon LAGRANGE
i+m (~) ~=H' h~
Fehlerabsch~itzungen und Extrapolation mit rationalen Funktionen 415
Die Formeln (3) und (5) gelten fiir beliebige hi, soferne nur Yi=/'Y ist. Fiir all- gemeinere Sequenzen yi ist eine dem Nevilleschen Algorithmus entsprechende, einfache Interpolationsformel nicht bekannt.
Im Falle allgemeiner Sequenzen yj .lassen sich jedoch die Operatoren A~ fiir die speziellen Schrittweiten hj=hoY, 0
416 ROLAND BULIRSCH und Jos~v STOER :
Satz t gibt fiber das Konvergenzverhalten dieser Extrapolationsverfahren Auskunff. Besitzt T(h) eine asymptotische Entwicklung (1) von h6chstens k Gliedern (k= oo ist m6glich), so konvergiert der Fehler in der n-ten Spalte yon (4) fiir n - -k+l wie h~ k+l. AuBerdem sieht man, dab die Absch/itzungen ftir den Extrapolationsfehler I~ Cm-i)- T(0)[, 1"= k, k + t . . . . . m, dieselbe obere Schranke liefern. Das legt die Vermutung nahe, dab nur die Berechnung der ersten k + l Spalten von (4) sinnvoll ist. Dies wird durch die Erfahrung best/itigt. Existiert ferner ftir jedes k das Restglied Rk+ l(h) in (t) und bleiben die Rk+l(h ) gleichm~i6ig beschriinkt, so zeigt Satz t, dab
ra ?m+t --1 ]T~ ) - T(0)I fiir m-+oo superlinear, n/imlich wie h, b~ -x~j gegen 0 strebt.
Zum Beweise yon Satz t werden zwei Hilfss~itze vorausgeschickt. Zun~ichst folgt fiir das Extrapolationsverfahren (3), das auf dem Neville-Algorithmus be- ruht, das
und hi4-1 ~b dann ist [iir k O, -[,j < 1,
~ I c~~ hl ~+ "~
FehlerabschAtzungen und Extrapolation mit rationalen Funktionen 41 7
Analog erh~ilt man
2 2 i=m-~ i=o h~-k+i-x l h~_a h~-k+ i h~-k+ i
ml m- -k - -1
U h~ x h~_~+~ #=m--k+~+l t /z~O 1
< h~m_~... ~ ~"~'~ k) < h2_~.., h~ C';(bV) : - t~mt- , l~ , , = 9
Die erste Behauptung ergibt sich jetzt mit
-C1 ( b~/, k) =C l ~,, k) 21- -,,C1 (b'V, k) ~ C 1 (bY).
Weiter gilt mit Hilfe der Produktdarstellung von sin z
X 1
~ __ (J+~ 't~[ j+~ , j+2 , , 1 ' \ rn~/ I " ] 2 " ~-+l - - '
und allgemein /Sir beliebiges k ~0 ganz,
4t8 ROLAND BULIRSCH und JosEF STOER:
Ftir q_>89 gelten die ftir diese Zwecke ausreichenden N~iherungsformeln
[v~3 (0' q)va4 (0' q)]-89 ~ .]/--~n~- exp ( - 2 ~ 8~q)
(2./,. l g--'nq 61nq
Beweis. Man hat
2 I c~}l by,+, : h~*+, y, l c(o~ (by*+,)J = ao~*+' (by*+, + bY=)... (bY*+, + bY,) "~ (l--b~,) (1 --bYe,) i=o i=o "" '
(1 +bYk+,-rk+x) ... (1 +by=-r*+,) . (t +br*+,-y*) ... (1 -[-byk+,-r,) = X (1 --bY,)... (1 --bvm)
k
h~k+l b I'~-k) rk+l+X=oyj
< (1 -- q2) ... (1 -- q2~) 2 (! q_qU)... (t +q,(~-k-x)) (t +q2) ... (1 +q2k) ( ! -q~) . . . (1 _q2(~+~))
hgk+, b("-k)v~+,+~Xo ~j.
Das Weitere folgt ]etzt aus den Beziehungen
Oi (O,q)=2q 88 3, a2(O,q)=2qt~l ( l+q2")2( l - -q ~') und n=l n=l
0;(0, q)=02(0, q) Oa(O, q) Oa(O, q).
Am Rande sei vermerkt, dab fiir y i= (l +])y , (q= b y12)
lim ~,[c~}[ = l"O--q=)'"O-q~O ,,,-+~o ~= o I ( l+q2) . . .O+q 2t) [Oa(O' q) O4(O' q)]-~"
In diesem Fall sind die Absch~itzungen in Lemma 2 scharf. Speziell ftir q= 89 (s. z.B. Romberg-Verfahren) ist
[0 a (0, ~) ~4 (0, ~)]- ~ = t ,969. . . und
02 - ~ = 5,3 . . . " s (0, ~-) , (o, q) ]
letztere Absch~ttzung ftir C2(~, k, 0) ist scharf, falls mi tm auch m- -k und k nach e~ strebt.
Der Beweis yon Satz t ist jetzt einfach. Ftir das Extrapolationsverfahren (3) ergibt sich aus Lemma I wegen hi
Fehlerabsch/~tzungen und Extrapolation mit rationalen Funktionen 4t 9
(8) liefert jetzt k
(m--k)~'k+l+ l~ '/j IT& ~ T(o)l ~M~+lC(bV)h~ ~+'b t=1 ,
womit der Beweis yon Satz I erbracht ist, denn der Ubergang yon (h o, T~ ~ nach (h i, T~ I) kann als , ,Umnumer ierung" der h s interpretiert werden.
Als letztes last Satz t die Frage nach dem Konvergenzverhalten von T~ ~ often, wenn von T(h) nut bekannt ist, dab lira T(h : )= T(0). Es gilt hier
]--+co
Satz 2. Es sei lira T(hi) = T(O). Ist /iir alle f entweder Yi =]~ und hi+~ < b < 1 i~oo hf-i - - oo
(Rekursion#ormel (3)) oa : hi= h#, 0 < b < ~ . .d Z= b" kon~.'g~nt (R~k.rsio~s- /ormel (7)), so gilt lirnooT2)= T(0), i =t , 2 . . . . -
Beweis. Ftir den Fall y /=]y ist dieser Satz bekannt ([1], [7], [3]). Ftir den Rest t folgt aus G1. (2a) nach TOEPLITZ (S. [3], Satz t) l i rn T~0= .lim T(h:), falls
l-+OO i - t -m
1. .2 4~ =1,
i+m
2. 21~l - - - -
420 ROLAND BULIRSCH und JOSEF STOER:
Das nAchste Beispiel betrifft die Integration der Differentialgleichung y '= y, y (0 )= t, mit Hilfe der Runge-Kutta-Formeln. In diesem Fall lautet die Ent- v cmung (vg l . E4] )
T(h, x) = y (x) + h" za (x) + h 5 ~ (x) + . . . + R (x, h).
Zur Extrapolat ion wurde fiir Am=2 -m G1. (7) mit b=~, Y i=3 +] verwendet.
Tabelle 2
#o.)= T(~.,X) ~=o)
o t 2 3 4
2,7083 3333 33t 2,71734619139 2,71820993 9t9 2,7t82 7684440 2,71828150036
2,71828182 844 (e = 2, 7182 8182 846)
2,7t 79 4704 860 2,7182 7786025 2,71828181 110 2,71828182844
,
Die bisher untersuchten Verfahren beruhten auf der Interpolation durch polynomartige Ausdrticke. Es liegt nun nahe, im Fall Y i=JY statt dessen ra- tionale Funkt ionen zu benutzen und als Niiherungswert ftir T(0) den Weft
T~(~)= Tv(~ ) (0) derjenigen rationalen Funkt ion in h v
Pu+- , +"" (9) ~.(o (h~- ~)~(h) !O ~!Oh:" +p~)hm,
zu nehmen, ftir welche gilt
(tO) T(v~(h j )=T(h i ) , j= i , i+ l . . . . . i+#+v=i+m.
In [19] wurde gezeigt, wie man die Werte "(0 T~, ~ (8) an einer festen SteUe 8 aus den gegebenen Werten T(hi) rekursiv berechnen kann, ohne die rationale Funkt ion
(i) Tvt~ (h) selbst, d.h. ihre Koeffizienten, zu bestimmen. Setzt man in die in F12] angegebenen Fonneln ~ = 0 ein, so erhalt man
To(i) _ T(hi) , o - - hVTfl+ a) -- h v 7"(0 , ' / ' ,(1)_ i /.*--1,o ;~+p .U--l,0
~,u,o--
und ffir/z, v ~ t To!~+xx) To!'~)-x
= - - ( -
h~T.(i+x) tT(i) _ T.(i+I) ~ _ h ~ T..(i) iT(i+ 1) _ T(i+I) T.(~) = i /~,V--ll /*,t'--I .U--l,t'--l,' i+ l~+V .u,~'-- I ~, I * ,V- -1 I~- - l , v - - l ! I,~,n h~(T(i)_ l_T(i+l) , ~, (T(i+I)_T{i+I) ~
- /a - - l , y - - l / - - '~+/~+v ~-p ,v - -1 IJ--l,v--l]
FehlerabschAtzungen und Extrapolation mit rationalen Funktionen 42t
W/ihlt man fiir die Extrapolation folgende Sequenz der (/z, v) : (0, 0) --> (0, 1) -+ (t, t) ---> (l, 2) ---> .-- -+ ( i , i ) -+ ( i , i+ l ) -+ ( i+1 , i+1) -+ .--. (Z/ihlergrad ---- Nennergrad (-- 1)) und schreibt man T(0,v+~ ftir T(0_~,,,, (# +v=m), so reduzieren sich die obigen Formeln auf
-(0 = 0 1
T~o O = T(h~)
T(i) -- T(i+ ~) m-- m--1 "~-
(1 t) T(i+ ~) _ Tj0
~ 2 ~ ) J --t
Man erkennt die forlnale Ahnlichkeit mit dem Neville-A]gorithmus (3)- Ver- wendet man die Formeln (11), so muB man a11erdings prtifen, ob die auftretenden Nenner nieht zu klein sind. Insbesondere besteht diese Gefahr bei gr6Beren m, da T. (i+1) T. (i+~) I ~- i - - ,~-~ t bald sehr klein wird. Dieses Verhalten ist jedoch ungeffihr- lich, weil es nur die schnelle Konvergenz des Verfahrens best~ttigt und in der Regel nicht auf algebraischen Ausnahmesituationen beruht, wie sie in [12] aus- fiihrlich beschrieben sind.
Bei der Extrapolation durch rationale Funktionen 1/iBt sieh nun ebenfalls ein Ausdruck ffir den Fehler T~(~ - T(0) angeben. Ist n/imlich kG#, so gilt wegen (9) und (10)
p lO (h.~ : T(hi) Q~!~ (hi) f : i, i + 1, i + t z + v : i + m. !a~vk 11 ~ " . '1
Multipliziert man diese Gleichungen mit den Lagrange-Koeffizienten c~}/ (s. (5)) und summiert fiber ], so erhltlt man wegen (1) und (2)
i+m
~(o) = ~(o) 0.~(o) + .X f~.h~ ~+~ &+~ (hj) O.~'),(hj) oder falls Q~), (0) = q(~) 4 = 0
, i+m
(t2) T(0 -- T(0) = " x-~ ~(i) L(/~+l)y/~ 9 _ . , , ~ ,.,,,;,oj ..~+.(h;) 9~!~ (hi). o ~=i
Leider enth/ilt dieser Ausdruck ftir den Extrapolationsfehler noch das Nenner- polynom Q(~I (hi) , so dab er zun~tchst nur bedingten Wert hat. WAre allerdings Rk+l(h i )=const , ~=i . . . . . i+m, so erg~be sich ftir k=#
# T(i)-u, ~ -- T(0) : ( - t)*~ h~ """ h~+" q-~00 Rk+t' I z+v:m"
Gilt lediglich lira T(hi)---- T(0), so konvergiert TJ~)~ ftir # +v-+ oo gegen T(0) ~--~ oo
(vgl. Satz 2), falls ~ ~ b < tund
(2~i')" (hi) < Ci, fiir alle /z, v, ]. (13) O~.~,(o) =
Denn wie in (t2) folgt
i+m i+m T(i) : t -~"* Q~),(O) ~' dO.n(O [h.~ T(hi) ~,3~.T(hi)
~=i =i Q~),(o) 29*
422 ROLAND BULIRSCH und JOSEF STOER"
und mit den c~} (s. [8], Satz t) erftillen wegen (t3) auch die ~-~ die Bedingungen von TOEPL ITZ
i+m
~-. =
Fehlerabsch~tzungen und Extrapolation mit rationalen Funktionen 423
zun~ichst den SchluB nahezulegen, dab k = 0 den kleinsten Fehler bei der Extra-
polation liefert. Andererseits gehen jedoch in q~k, " wie sich zeigen l~tl3t, qT Ableitungen bis zur Ordnung/(2,+~) (x), r= max {k, m- k}, in komplizierter Weise ein. Man wird also in der Regel am gtinstigsten mit k=/gm/ (Z~ihlergrad ~
r 1
L~3
Nennergrad) arbeiten. Die numerische Erfahrung auf der Rechenanlage PERM (TH Mtinchen) hat diesen SchluB best~itigt.
In den folgenden Beispieten zur Trapezsummen-Extrapolation sind die Er- gebnisse bei unterschiedlicher Extrapolation gegeniibergestellt. Als Sehrittweiten-
folgewurde die bereits in [31 erw/ihnteFolge~ {L, L2 ,3 ,4 ,6 ,8 , t2L L L L L . . . . } gew~thlt (L L~tnge des Integrationsintervalls), die auch hier in bezug auf Arbeits- aufwand und Genauigkeit die besten Resultate liefert. In der folgenden Tabelle gibt m die Anzahl der Extrapolationsschritte, A,~ die Anzahl der bis zur Schritt- weite h~ zu berechnenden Funktionswerte und T(hm) die berechneten Trapez- summen an. Die weiteren Spalten enthalten in der Reihenfolge die extrapolierten Werte bei der Interpolation durch
a) Polynome,
b) rationale Funktionen (Z~ihlergrad = Nennergrad (--!)), c) reziproke Polynome.
(Es wurden zwei ,,sehlechte" und ein ,,gutes" Beispiel ausgew/ihlt l)
Tabelle 3
m
e
1
2
f sin x dx 0
- - f xsin3xdx 6
Am T(hm)
1,00585652083 1,00330707022 1,0014 7397628 1,00082994518 1,00036913108
0,97704861 6664 o,9871t58oo974 0,994281888297 0,996785171887 0,9985716979o7
0,2741 5567 7808 0,249811161476 0,23400650 8335 0,228760775525 0,2251 00165925 0,223835594426 0,22293754 7635
~!o
1,0000 0281 498 1,0000 0014667 1,0000 0000488 1,0000 0000010 1,0000 0000000
0,0999984959951o~; t,00000005825 0,09999999999210, 1,0000 0000000
0,1882 34799393 0,223744283778 0,222t 93662180 0,222222514134 0,222222220926 0,222222222186 0,222222222213
T'(~ 2 [m121, m- - Eml ]
t,00000100 964 1,00000003 995 1,00000000080 1,00000000000
0,0999 9957021010, t,00000000656 1,0000 0000000
0,2226 5227 5416 0,22222326 8686 0,22222221 7815 0,222222222178 0,2222 22222223
T(O) o ,m
1,00000356225 1,000000t7065 !,00000000541 1,00000000010 !,00000000000
i 0,9994 8556 0817 1,00000816914 0,9999 9994 2886 t ,0000 0000 021 t ,0000 0000 000
0,2223 4864 0655 0,2222 2334 67O8 0,2222 2222 5452 0,2222 2222 2186 0, 2222 2222 2223
Man sieht, dab die Werte T, tO}~l,m_E,~/21 jeweils am schnellsten konvergieren. Ein weiteres interessantes Beispiel liefert die Integration der Differential-
gleichung y'=/(x, y) nach dem Euler-Verfahren y,+l=y,+h/(x, , y,). Das
424 ROLAND BULIRSCH und JosEF STOER:
Fehlerverhalten ist bier gegeben durch
(t5) T(h, x)=y(x) + hzl(x) + h***(x) + Mz3(x) +....
Als spezielles Beispiel wurde wieder y'= y, y (0) = 1 gew~hlt, vgl. dazu auch [1], w (In diesem Fall ist T(h, x)=e*--h2eX+h~(3 + x-~)e~+ ...).
Mit der Schrittweite h,~=2 -m ergibt sich
Tabel le 4
T(hm, 1 ) T(mO~O T[(~/21, m -- [m/~]
2,56578451393 2,63792849736 2,67699012937 2,6973 4495258
2,71563200012
2,7t38 7899486 2,718o2983460 2,71827434380 2 ,7 i82817t 514
2 ,7 i828182855
2,7t81 55 i6 t73 2,71828078589 2,71828181 79i 2,71828182855
Man erhitlt also mit der einfachen Eulerschen Methode for m = 6 dutch Bildung yon T~~ praktisch dieselbe Genauigkeit wie mit den Runge-Kutta-Formeln mit anschlieBender Extrapolation (vgl. das Beispiel am Ende von w 3). Die Anzahl der zu berechnenden Funktionswerte /(x~, Yn) betriigt in diesen beiden F~tllen 121 bzw. t20. Das Verh~iltnis verschiebt sich sogar noch zugunsten der Euler- Integration, falls die Schrittweitenfolge ~ benutzt wird. In diesem Fall erhiilt man durch Berechnung yon nur 30 Funktionswerten das Ergebnis T~,~ 2,71828182652. H6here Genauigkeiten k6nnen hier mit ~ im Gegensatz zur Trapezsummen-Extrapolation leider nicht erzielt werden, da wegen der lediglich nach Potenzen yon h fortschreitenden Entwicklung (15) dutch die nahe bei 1
9 hi+l die liegenden Quotienten~--. " Rundungsfehler stark ins Gewicht fallen. OG
Das folgende AzGoz-Programm zur Berechnung von f/(x)dx benutzt zur UG
Trapezsummen-Extrapolation rationale Funktionen (Z~thlergrad ~ Nennergrad,
s. Formel (1t)). Verwendet wird die Folge ~={ L'L2'3'4'6'8L L L L . . . . }, L= OG--UG. Das Programm ist optimal konstruiert, bereits berechnete [(n hi) werden zur Ermittlung neuer To C'~) wieder verwendet. Bezeichnet A m die Anzahl der [(n hi) zur Berechnung von To ~~ . . . . . To (ml, so ist
2 -+l , m gerade Am = t + | m+l m-1
~2 ~- +2 ~- , m ungerade
m_>_2
(vgl. Tabelle 1 in E3]).
Fehlerabschi~tzungen und Extrapolation mit rationalen Funktionen 425
Das Programm ermittelt zu vorgegebenem ord die T~(0_ i ffir
m=0, t, . . . , ord
i = m, m -- 1 . . . . . r (m) mit r (m) = max {0, m -- 7} Oa
und liefert als NAherung ffir f / (x )dx den Wert To(~'~o~)al. va
(Tocie Erfahrung zeigt, dab es keinen Sinn hat, fiber mehr als 8 Stfitzpunkte (hi, 1"1) zu extrapolieren.)
Voraussetzung ffir die schnelle Konvergenz der T~_i ist allerdings die Existenz OG
einer Entwicklung (1) ffir T(h) mit 7i-----2j und k>=7 (trifft zu, falls f [/~ze)(x)] dx existiert), u a
reol procedure Trapezsummenextrapolation (ug, og, ord, eps) procedure : ([) ; comment Die Prozedur Trapezsummenextrapolation liefert den N~herungswert
ffir das Integral der Funktion [ (x) zwischen der unteren Grenze ug und der oberen Grenze og, den man durch ein Extrapolationsver- fahren mit rationalen Funktionen vonder Ordnung oral (>--_ 2) erh~ilt. Die Extrapolation wird vor ihrem natfirlichen Ende abgebrochen, wenn sich zwei aufeinanderfolgende Niiherungswerte t [i] und t [i + 1] ffir das Integral um weniger als eps abs (t [i + 1]) unterscheiden;
volue ug, og, ord, eps ; real ug, og, eps ; integer ord; real procedure [; begin real h, e, tO, t2a, t2, tn, tS, ha, hg;
integer m, nn, i; integer arroy n [0 : ord I ; arroy t [0 : 71 ; boolean bo ;
procedure extr (m) ; value m, integer m; begin real v, d, u, by, hu;
integer i, mr; v:---0; u:--rE0]; h:----rE0] :----tn; if m> 7 then mr:= 7 else mr:=m; for i:----t step I until mrdo begin d :=n[ml /n [m- - i ] ; d :---d
hv :=h- -v ; hu : :h - -u ; if hv :~0 then
begin h := h+hu/ (dx ( l - -hu/hv) - -1) ; v :=u; u :=t[ i ] ; t[i] :=h;
end else go to ende;
end; end ex~r ;
426 ROLAND BULIRSCH und JOSEF STOER :
n[O]:=1; n[I] :----- 2; n [2] :=3; fo r i : ---- 3 s tep I unt i l ord do n [i] : = n [i - - 2] 2;
e :-~ og- -ug ; bo := true; tO :----- ([(og) +/ (ug) ) /2 ; [ [0] : = tO e; tea : = te : = [ (ug + e/2) ; tn : = (~0 + t2) e/2 ; extr (t) ; t3 : = / (ug + e / 3) + / (og - e / 3) ; tn : = (tO + tB) e / 3 ; extr (2) ; ha : = h;
for m := 3 step I until ord do begin nn :---- n [m]; hg :---- e/nn;
if bo then beg in fo r i :---- 1 s tep 2 unt i l nn do
t2:---- t 2 + [ (ug + i x hg) ; tn : = (~2 + tO) hg;
end else
begin for i := t step 6 until nn, i:---- 5 step 6 until nn do tB := tB + /(ug + i xhg) ;
tn : ---- (t3 + t2a + tO) hg; t2a : ~- t2; end;
extr (m) ; bo : = 7 bo ; i f abs (ha - - h) < abs (h) x eps then go to ende; ha :----h;
end; ende : ~rapezsummenextrapolation : ---- h; end trapezsummenextrapolat ion
Li teratur
[1] BAUER, F. L., H. RUTISHAUSER, and E. STIEFEL" New Aspects in Numerical Quadrature. Proc. of Symposia in Applied Mathematics 15, Am. Math. Soc. (1963).
[2] BOLTON, H. C., and H . I . Scoll~S: Eigenvalues of Differential Equations by Finite-Difference Methods. Proc. of the Cambridge Phil. Soc. 52, 2t5--229 (1956).
[3] BULIRSCH, R.: Bemerkungen zur Romberg-Integration. Num. Math. 6, 6--16 (t964).
[4] GRAGG, W. : Repeated Extrapolation to the Limit in the Numerical Solution of Ordinary Differential Equations. Thesis UCLA (1963).
[5] HARTREE, D. R.: Numerical Analysis. Oxford: At the Clarendon Press 1958. [6] KUNTZMANN, J.: M6thodes num6riques, interpolation-dgriv6es. Paris: Dunod
1959. [7] LAURE~T, P.- J . : Un th6or~me de convergence pour le proc6d~ d'extrapolation
de RICHARDSON. Comptes Rendus de l'Aead6mie des Sciences (Paris) 256, 1435--1437 (t963).
[8] RICHARDSOn, C., and J. GAUNT: The Deferred Approach to the Limit. Trans. Roy. Soe. Lond. 226, 300--361 (1927).
[9] ROMBERG, W. : Vereinfachte numerische Integration. Det. Kong. Norske Viden- skabers Selskab Forhandlinger 28, Nr. 7, Trondheim t 955.
[10] RUXlSHAUSER, H. : Ausdehnung des Rombergsehen Prinzips. Num. Math. 5, 48--54 (1963).
Fehlerabsch~tzungen und Extrapolation mit rationalen Funktionen 427
[11] STETTER, H.-J.: Asymptotic Expansions for the Error of Discretisation Algo- rithms for Non-linear Functional Equations. Erscheint in Num. Math.
[12] STOER, J. : Uber zwei Algorithmen zur Interpolation mit rationalen Funktionen. Num. Math. 3, 285--304 (1961).
[13] COREY: The American math. Monthly 13 (1906). [143 SALVADORI, M. G., and M. L. BARON : Numerical Methods in Engineering. Engle-
wood Cliffs, N.J.: Prentice-Hall 1952. [15] LVNESS, J. N., and B. J. J. McHuoH: Integration over multidimensional Hyper-
cubes, I. A progressive Procedure. The Computer J. 6, 264-270 (1963). [16] LAURENT, P.-J. : Formules de quadrature approch6e sur domaines rectangulaires
convergentes pour route fonction int6grable Riemann. Comptes Rendus de l'Acad6mie des Sciences (Paris) 2S8, 798--801 (1964).
[17] NAVOT, I.: The Euler-Maclaurin Functional for Functions with a Quasi-Step Discontinuity. Math. Comp. 17, 337--345 (1963).
Mathematisches Institut der Technischen Hochschule 8 Miinchen 2, Arcisstr. 21
(Eingegangen am 26. Mtirz 1964)