ROCZNIKI POLSKIEGO TOWARZYSTWA MATEMATYCZNEGO Seria III: MATEMATYKA STOSOWANA XXIX (1987)
J
a n u s zC
h o j n a c k i, A
n d r z e jK
ie ł b a s iń s k i(Warszawa)
Błędy zaokrągleń w algorytmie Romberga
(Praca wpłynęła do Redakcji 1985.07.04)
0. Wprowadzenie
0.0. Ekstrapolacyjny schemat Neville’a (Richardsona) stosowany dla obli- czania kwadratur w metodzie Romberga, p. [2], [5], budzi niekiedy pewne obawy. Ekstrapolacja nie cieszy się bowiem dobrą opinią u numeryków.
Wiadomo, że może powodować istotną utratę dokładności.
W przypadku metody Romberga ekstrapolacja jednak szybko „wygasa”, gdy przechodzimy do obliczania kwadratur wyższych rzędów. Pozwala to żywić nadzieję, że utrata dokładności będzie niewielka. Eksperymenty nume- ryczne uzasadniają zresztą w pełni ten optymizm.
Pozostają jednak nadal pytania:
— jak można scharakteryzować dokładność kwadratur, obliczanych w metodzie Romberga?
— jaką część ewentualnie utraconej dokładności należy przypisać stoso- waniu ekstrapolacji, a jaką prostej kumulacji błędów zaokrągleń?
— od jakich etapów obliczeń zależy w głównej mierze dokładność otrzy- manych kwadratur?
— czy można kumulację błędów uniezależnić od wymiaru zadania?
W poniższej pracy próbujemy odpowiedzieć na te pytania.
0.1. Przypomnijmy tu, że każda z kwadratur Romberga jest wyrażeniem postaci
n n = i B j ( x k),
k= 1
gdzie f oznacza funkcję podcałkową, a układ węzłów i wag kwadratury T, {**, Bk}k=l, definiuje ją jednoznacznie.
W naszej analizie za dane początkowe algorytmu Romberga przyjmujemy obliczone wartości {fk}k= i funkcji f w węzłach X k.
Błąd reprezentacji liczby f k w arytmetyce numerycznej, por. [6] lub [2], nie
przekracza wielkości v |/k|, gdzie v oznacza tzw. względną precyzję arytmetyki
(w typowych arytmetykach współczesnych zachodzi 10- 14< v ^ 10"6). Za-
26 J. Chojn acki, A. Kieł bas iń ski
tem poziom błędów reprezentacji danych w kwadraturze T [ /] daje się oszaco- wać wielkością
P := v £ \Bk-fk\.
k = 1
W pracy pokazujemy, że wpływ błędów zaokrągleń (występujących w algorytmie Romberga) na obliczoną kwadraturę daje się oszacować wielkoś- cią f P, gdzie f zależy od liczby N węzłów kwadratury, f = f (N).
Rozróżnienie w analizie zjawisk: prostej kumulacji błędów oraz ekstrapola- cyjnego powiększenia błędów, prowadzi do przedstawienia mnożnika f w postaci sumy:
T =
XK+
XE,
przy czym „składnik ekstrapolacyjny” xE jest zawsze rzędu jedności, nato- miast „składnik kumulacyjny”, xK = t a :(N) jest rzędu N, gdy w algorytmie stosujemy proste sumowanie (wyraz za wyrazem), a rzędu log 2 N, gdy stosujemy bardziej wyrafinowany algorytm sumowania.
Przedstawiamy wreszcie wariant algorytmu Romberga, w którym f nie zależy od N.
0. 2. A oto zawartość dalszych rozdziałów pracy. W rozdziale 1 przypomi- namy prosty algorytm metody Romberga i jego własności istotne dla dalszej analizy. Przeprowadzamy także wstępne badanie procesu powstawania i przenoszenia błędów zaokrągleń w tym algorytmie.
W rozdziale 2 definiujemy modelową arytmetykę numeryczną i badamy błędy zaokrągleń w prostym algorytmie Romberga. Uzyskujemy tu podsta- wowe twierdzenie, charakteryzujące dokładność obliczonych kwadratur.
W rozdziale 3 analizujemy stabilność badanego algorytmu, podajemy oszacowania błędu względnego obliczonych kwadratur. Rozważamy modyfi- kacje oszacowań błędów przy założeniu regularności funkcji podcałkowej.
Dyskutujemy konsekwencje stosowania różnych algorytmów sumowania oraz różnych arytmetyk numerycznych.
W rozdziale 4 rozważamy inne\varianty algorytmu Romberga. Wskazuje- my wariant, w którym kumulacja błędów nie zależy od liczby węzłów kwadratur.
W rozdziale 5 przedstawiamy informację o własnościach i implementacji ważniejszych algorytmów sumowania.
1. Algorytm Romberga, wstępna analiza jego dokładności
1.0. Kwadratury Romberga. Rozważamy funkcję / określoną na przedziale
[u, 6] oraz kwadratury Romberga {7}(,)} tej funkcji, generowane przez poło-
wienie kroku dyskretyzacji (por. [2] rozdz. 4, [5] rozdz. 3).
Błędy zaokrągleń w algorytmie Romberga
27 Kwadratury te tworzą tablicę trójkątną
t u )
73'01 T
q1) Tll) Ttf) 1o 7^(0 77(0
w = T ił'm = b- ^ Y . " /(*■?), Pierwsza kolumna tej tablicy zawiera kwadratury trapezów
(
1
.2
)z węzłami (1.3)
k = 0
... b — a
= a-\~k— —, k — 0, . . . , 2‘. 2 ’
Symbol sumowania z dwoma primami wyjaśniony jest wzorem
(1.4) Z rt @
qa k ~ ■ y + fli + ••• + a N - l + — . k = o
Zauważamy, że kolejna kwadratura trapezów może być obliczana z poprzed- niej przez dołączenie odpowiedniej krotności sumy wartości funkcji podcał- kowej na nowych węzłach:
(1-5) 7y> = l 73‘- » + ^ J£ /(x S -,).
Aby to uzasadnić, wystarczy skorzystać z (1.2) oraz uwzględnić zależność (por. (1.3))
(1.6) xS = x r i), 1 = 0 , l , . . . ^ ' " 1.
Kwadratury j-tej kolumny tablicy (1.1) wyrażają się przez kwadratury kolu- mny (7 —l)-ej ekstrapolacyjną zależnością:
4jrli) _ Tli - 1)
(1-7) 7 f = ■ 1
Wychodząc z tej zależności, można wykazać, że
( 1 . 8 )
L k = 0
przy czym wagi f}jk spełniają związki rekurencyjne
'>-9) = ( * P j- 1,* - 2fij-..V2)/(4J- O.
28 J. Chojn acki, A. K iełb as iń sk i
gdzie
( 1 . 10 )
Po,k ~1>
P j - i,k/2Pj-i.i, gdy k = 21, 0, gdy k = 21 — 1
(przy / — naturalnym). Przez prostą indukcję otrzymujemy następnie (por.
[ 1]):
( 1 . 11 )
P j,2 i
-i =
P j,i = ń
P j.2 =5 ( 1 + 2/4')•
PM, 0.484 < Pj'2 ^ Pjtk < PjA < 1.452.
Zależności (1.5) i (1.7) stanowią podstawę konstrukcji następującego algorytmu.
1.1. Algorytm Romberga. Algorytm oblicza m +1 pierwszych wierszy tab- licy trójkątnej {Sf}i=0....my=0,...,o takiej, że
(
1
.12
)Danymi początkowymi są ciągi {f&0), f i 0)}, {/2(,?-i}‘ = J....gdzie
(1.13) / k(i> = (wartość funkcji f obliczona w węźle b — a a + k -— r } .
Algorytm 1
(1.14) sj>0>: = ( / r + / i (0))/2;
for i : = 1 to m do begin
(1.15) SN(i):= Z f2 ł-u (1.16) $<*>:= +
for j:= 1 to i do
(1.17) S f: = (4^'^S’yi i - S ^ / t y ^ - l )
end; □
U w ag i, (i) Obliczanie elementów tablicy (S)0} zamiast {77°} motywowa- ne jest pewnym zmniejszeniem liczby działań oraz zaokrągleń, jakie przez to osiągamy. Spośród ~ m 2/ 2 kwadratur {7}(,)} potrzebujemy zazwyczaj tylko jednej, którą możemy obliczyć z odpowiedniego S{j\ stosując instrukcję (1.12).
Zaliczmy tę instrukcję do algorytmu 1, tak jakby była wykonywana dla wszystkich i = 0, 1,..., m, j = 0,..., i.
(ii) W realnych obliczeniach na ogół nie musimy posługiwać się tablicą o
m2/2 elementach. Ponieważ obliczamy elementy tablicy {S)ł)} wierszami, więc
na ogół wystarcza przechowywanie elementu S f w zmiennej S [7] jednowy-
miarowej tablicy S[0:m ] (jak długo ta wartość nie będzie zastąpiona S{j +i)).
Błędy zaokrągleń w algorytmie Romberga
29 1.2. Wstępna analiza powstawania i przenoszenia błędów w algorytmie 1.
Niech S f oznacza wartość obliczoną, a Ś f wartość dokładną (która powsta- łaby w algorytmie, gdyby był realizowany bez błędów zaokrągleń). Interesują nas błędy A f — S f —Ś f. Wszystkie liczby S f są przybliżeniami tej samej
b
wielkości j f (x) dx/(b — a), więc sensowne jest porównywanie błędów {Af}
między sobą.
aPonieważ obawy budzi „etap ekstrapolacyjny” obliczeń, więc prześledzi- my jak powstają i przenoszą się błędy przy wykonaniu instrukcji (1.17).
Zachodzi równość
(1.18) S f = Ś f + A f = 4^(Ą»1+ 4 i 1)-(Ą L -11,+ 4 - i 1)) 4j — 1 + tj,(0
gdzie y f oznacza błąd, wytworzony zaokrągleniami w instrukcji (1.17).
Otrzymujemy stąd zależność
(1.19) ... 4 j A(j l l - A (i l ł ) ...
J 4j - 1 Y} ’
w której pierwszy składnik nazwiemy błędem przeniesionym, a drugi — błędem wytworzonym. Dla błędu przeniesionego uzyskujemy oszacowanie:
4 ' A f - i - A p P
4j-
* r -t+
A d) 1)
A j - 1 A j - I
4J — 1
<< ( 1+ ^7T Y )m ax( l4 - il» 1 4 - ^D*
A więc już samo przeniesienie błędów z ( j— l)-ej kolumny na j-tą może powiększyć poziom błędu z mnożnikiem ^1 Przyczyną tego powięk- szenia poziomu błędu jest oczywiście stosowana ekstrapolacja. Zapisując (1.19) w postaci
(1.20) A f ^ A f . i + a ł f + y f, możemy składnik
(1.21) ^ > = 3 - 1 , AJZ1.
nazwać przyrostem ekstrapolacyjnym błędu, a składnik y f — przyrostem kumulacyjnym.
Zauważamy przy okazji, że przy przechodzeniu do dalszych kolumn (tzn.
przy rosnącym j ) znaczenie przyrostu ekstrapolacyjnego szybko maleje („eks-
trapolacja szybko wygasa”), zatem już ta wstępna, pobieżna analiza pozwala
30 J. Chojn acki, A. Kieł bas iń ski
uciszyć obawy, związane ze stosowaniem ekstrapolacji. Dopiero jednak bar- dziej szczegółowa analiza z rozdz. 2 i 3 wykazuje, jak nieznaczne są straty dokładności wywołane ekstrapolacją w tym procesie.
2. Błędy zaokrągleń w prostym algorytmie Romberga
2.0. Arytmetyka zmiennopozycyjna (fl). Rozważamy tu rzeczywistą binar- ną arytmetykę zmiennopozycyjną z t-cyfrową mantysą, z poprawnym za- okrąglaniem wyników działań arytmetycznych, por. [ 2] rozdz. 1, [6] rozdz. 1.
Będziemy ją nazywali arytmetyką fl.
Zakładamy, że w badanym procesie nie występują nadmiary ani niedo- miary. Zatem dla dowolnego działania <>e{ + , —, *, /} i dla argumentów x, y efl obliczonym wynikiem jest
(2.1) z = fl(x <>y) = (x O y )(l+ e), kl < v.
Wielkość £ = e(x, y, O) jest tutaj błędem względnym zaokrąglenia wyniku rozważanego działania, a v:= 2~t jest względną precyzją arytmetyki fl.
Zakładamy również, że występujące w algorytmie liczby całkowite (np. 2‘, 4‘, 4* — 1) są obliczane dokładnie i mogą występować w działaniach arytmety- ki fl jako argumenty (na równi z liczbami tej arytmetyki), zgodnie z regułą
(
2
.1
) .Przyjmujemy ponadto, że mnożenie lub dzielenie liczby arytmetyki fl przez całkowitą potęgę liczby 2 jest wykonywane bez zaokrągleń, np.
(2.2) fl (x/2i) = x/2\ fl (4 j*y) = 4 j y.
2.1. Linearyzacja oszacowań błędów. Liczba v = 2~' jest we współcześnie stosowanych arytmetykach bardzo mała (10-14 ^ v ^ 10“ 6). Pozwala to na (przybliżoną) linearyzację niektórych wyrażeń, zależnych od v.
Jeśli dla niezależnych wielkości e i <p słuszne są oszacowania
(2.3) |e| ^ ve, \(p\ ^ vę,
z niezbyt wielkimi e i ę, to dla wielkości y, zdefiniowanej wzorem (2.4) 1 + 7 = (1 + e)(l + <p),
będziemy stosowali zapis równości przybliżonej
(2.5) y = e + (p
oraz oszacowanie przybliżone
(2.6) lvl < vy, y = £ + q>.
Taka linearyzacja jest oczywiście dopuszczalna tylko w przypadku, gdy zachodzi v2E(p v(e + (p). Nierówność taką nazywamy warunkiem linearyzacji wielkości y, określonej w (2.4). Warunek ten zapisujemy zwykle w postaci nieco ostrzejszej
(2.7) vy 1.
Błędy zaokrągleń w algorytmie Romberga
31 U w aga. Postulowaną niezależność t: i <p należy rozumieć następująco:
istotne części r. i ip nie zależą od tych samych wielkości.
Skrajnym przykładem zależności jest: e = — (p, e = ćp. W tym przypadku równość przybliżona (2.5) przechodzi w ,, - e 2 = 0”, zaś nierówność (2.6) w
„£2 ^ v • 2e ”, co jest mylące bądź nierealistyczne.
2.2. Dane początkowe algorytmu Romberga. Zgodnie z opisem algorytmu 1 danymi początkowymi są w nim liczby / k(,), tzn. obliczone wartości funkcji / w węzłach x[(). Obliczenie to musi być wykonane jakimś algorytmem, wybranym przez użytkownika algorytmu Romberga. Wytworzone przy obli- czaniu f ( x (^) błędy,
(2.8) dfP
wpływają oczywiście na dokładność obliczonej kwadratury, ale nie są przed- miotem naszych rozważań (gdyż nie są błędami zaokrągleń w algorytmie Romberga). Odwołamy się do nich jednak w rozdz. 3 przy okazji interpreta- cji wyników przeprowadzonej analizy błędów zaokrągleń.
2.3. Algorytm sumowania. W algorytmie 1 nie jest określony sposób realizacji instrukcji (1.15):
SN''':=2£ A Z= 1
W rozdziale 5 przedstawiamy kilka ważniejszych algorytmów sumowania, charakterystyki ich kosztu i dokładności. Charakterystykę dokładności wy- korzystamy w analizie błędów w następnych punktach tego rozdziału.
2.4* Analiza błędów zaokrągleń algorytmu 1
A. Obliczanie elementów pierwszej kolumny: Realizacja instrukcji (1.14) daje obliczoną wartość
s f = ( / r + / r m + f ? o ) / 2 ,, ieoi < v.
Aby uzyskać bardziej jednolity układ oszacowań błędów, popełnimy w tym miejscu drobne odstępstwo od przyjętych reguł zapisu, przyjmując równość
(2-9) Sg>,=C/3°>+/1(0,)/2,
w której f f 0) oznacza naprawdę f {0)( \ +
q0), i = 0, 1. Można to np. interpre- tować jako „włączenie” błędu f i 0) • g0 do błędu ófi0) dla i = 0, 1, por. (2.8).
Realizacja instrukcji (1.15) daje, uwzględniając (5.1)-(5.3), 2' “ 1
(2.10) SNm = Y. + W ’1 < >V»
1= 1 a instrukcji (1.16), por. (2.2),
(2.11) S f = (Sg“ ll/2+SN"'l/2'')»(l + e(), |ft| ^ v.
32 J. Chojn acki, A. Kieł basi ńsk i
Przyjmując założenie indukcyjne
2
i_
1S (o 1 ) = ^ T Z f k X)( l + e o,k X))» \ 4 , k 1]\ ^ v - e g 1), k= 0
otrzymujemy z (2.10) i (2.11)
( 2 . 12 ) S0 = 4 Z ” /*° (1 + 4 lk), l£o!fcl < vejj\
^ fc= o gdzie, por. (2.3)-(2.5)
W o ł k =f X k ] + Q h
(2.13) j£(ó , 1), gdy k = 21,
& > / , (o gdy , , k = 21— 1 J tót oraz, z uwzględnieniem (2.9)
|Śo0) — 0, £<‘> = ^ + 1, t f 0 = max(ejj-1), ę (i)).
Stąd oraz z tablicy 5.1 i (5.4) wyprowadzamy zależność %l) od i dla różnych algorytmów sumowania:
(2.14)
2.1. Tablica oszacowań (,.<>)) <fc0 )
Algorytm sumowania j?0
wzw 2‘~ 1 — 1 2 ' ~ 1
pdw, gm, gm2, wzw f|2 i — 1 i
B. Obliczanie dalszych kolumn: {S^}, j > 0. Realizacja „ekstrapolacyjnej’
instrukcji (1.17) daje obliczoną wartość elementu j-tej kolumny, por. (2.2), (2.15)
Ą j C(i) _ C ( « - 1 )
S P - ^ ( 1+'?“')■ W < v 2.
Przy założeniu indukcyjnym (dla q = i — 1, i):
1 2 «
(2.16) spi, = - Z +
Z k = O
otrzymujemy z (2.15)
(2.17) 1
^ k = O
gdzie, por. (1.9)-(1.10),
(2.18) + +
Błędy zaokrągleń w algorytmie Romberga
33
(2.19) 1,2/'
e(« - D lk = 2l,
k = 2/ - 1.
Możemy więc napisać
( 2 . 20 )
gdzie
(2.21) Ą!i = I {&, <3° = I
p = l p = 1
Odpowiednio, dla wielkości £^° z (2.17), możemy podać wzór
(2.22) ę = ® + BJ + 2j,
jeśli zachodzi
(2.23) I ® « v6j.
W następnym punkcie zajmiemy się wyznaczeniem oszacowania typu (2.23).
Tu zauważmy jeszcze, że (2.18) przedstawia rozbicie zaburzenia eJJ na:
— zaburzenie £j‘ł lk — przeniesione z poprzedniej kolumny,
— przyrost Cj!k zaburzenia w wyniku ekstrapolacji,
— przyrost rjf zaburzenia o nowe błędy zaokrągleń.
Zatem Q(£k w (2.20) stanowi „składnik ekstrapolacyjny” zaburzenia.
C. Oszacowanie składnika ekstrapolacyjnego zaburzeń. Z (1.9) i (1.11) wynika nierówność
(2.24)
P p - i , / < P p - 1,1(4P— l)Pp,2i ^ (4P— l )/?p>2 Zatem z (2.19) otrzymujemy
3 4p + 2 ’
(2.25) gdzie (2.26)
ICp+ 1,2/1
4P i + 2 l^p.^/l
,.(/) _ P(i- 1)
*'P, 2/ ‘*p,/
Ze względu na (2.20) zachodzi równość
<2-27)
a z (2.13) wynika
(2.28) ^0,21 ^ Qi, l^o!2/l ^ v-
- Matematyka Stosowana
34 J. Chojnacki, A. K iełb asi ńsk i
Wprowadźmy liczby {<5P}' = 0> {CP+i}i=o wzorami
(<5o = C p + i = Ą p + i _j_2 ^ p ’
(2'29) ■ p p+1 = S p + 4 + 2fp ł l .
Z (2.25), (2.27), (2.28) wynika, że zachodzą nierówności (2.30) \df'2,\ vSp, |C®t| <
Eliminując 5P, otrzymujemy z (2.29) zależność rekurencyjną dla {£p}:
(2.31) Ć! = 1, lP+1 = [(4p-t- 14)£p + 24]/(4p+1 + 2).
A oto tabelka kilku pierwszych wyrazów ciągu {£p} oraz jego sum częścio- wych
(2.32) 8j = t l t .
P= 1
2.2 Tablica oszacowań {#,}
j Ó
0 0.000
1 1.000 1.000
2 2.333 3.333
3 1.424 4.802
4 0.524 5.281
5 0.161 5.442
6 0.047 5.489
7 0.013 5.502
00
0 < 5.507
Można wykazać, że oszacowania 01? 02 są realistyczne, natomiast dla j > 2 w nierówności (2.23) prawa strona jest zawsze istotnie większa niż lewa.
Rezygnujemy tu jednak z bardziej złożonych rozważań, pozwalających po- prawić oszacowanie 03, itd., gdyż uzyskane grube oszacowanie (2.33) |0J?2|| ^ v*ą < v*5.5,
w zupełności wystarcza dla naszych dalszych celów.
D. Błędy w obliczonych kwadraturach. Obliczenie kwadratury 7}(i) ze wzoru (1.12) wiąże się z ewentualnymi dwoma dalszymi zaokrągleniami. Dla obliczonej wartości otrzymujemy równość:
(2.34) = ( b - a ) - S f- { \+ x f) , \ x f \ ś v-2.
Uwzględniając (2.17), (2.20), (2.13), otrzymujemy więc następujące
Błędy zaokrągleń w algorytmie Romberga
35 T
w ie r d z e n ie1. (i) Algorytm 1 (wraz z instrukcją (1.12)), realizowany w arytmetyce fl daje obliczone wartości f j l), O ^ j ^ i < m, kwadratur Romberga, przy czym istnieją zaburzenia (ijJJ, k = 0,..., 2‘, ta/cic, że:
(2.35) (2.36) (2.37) Zatem
j y - V z k=0 itjii < »•*?..
t JI = Cxli>+ + Ce<+ + «1°].
= + !>]] + & + 3]- 0381 I (Ul < „ }t2‘ +2/+ 7.5) dla sumowania wzw,
1(1 + 27 + 7.5) dla sumowania pozostałymi algorytmami (rozdz. 5).
(ii) Oszacowania poszczególnych składników z (2.36) otrzymujemy z na- stępujących zależności (lub tablic):
(2.39) .(i-1)
’0,lj
W(0
* = 21, = O, +
gdzie <p(,) określone są w tablicy 5.1, a „składniki ekstrapolacyjne” Ofk zaburzeń spełniają zależności
(2.40) 0j?2/-i = 0, |0S?2il^v0„
gdzie Oj są określone w tablicy 2.2. □
3. Numeryczna stabilność algorytmu Romberga
3.0. Stabilność algorytmu. Z twierdzenia 1 wynika, że obliczona wartość T}(1) jest dokładną kwadraturą dla zaburzonych danych
(3-l) +
ze względnymi zaburzeniami Oznacza to 'numeryczną poprawność (ze względu na dane początkowe {/fc(,)}) algorytmu obliczania 7}(,) z charakterysty- ką kumulacji fj0:
(3.2) 2I_1 + 2j + 7.5
i + 2j + 7.5
(sumowanie „wzw”), (sumowanie „pdw” i inne).
Ponieważ zaburzenia xfk są odmienne dla każdej obliczonej kwadratury 7)°, więc algorytm obliczania całej tablicy {7](,)}i=0....mJ=0...., jest tylko numery- cznie prawie poprawny (a więc numerycznie stabilny, por. [2], rozdz. 1) z charakterystyką kumulacji K m:
(3.3) c m
_2m_1 + 2m + 7.5 3m + 7.5
(sumowanie „wzw”),
(sumowanie „pdw” i inne).
36 J. Ch ojnacki, A. K iełb asi ńsk i
Oszacowania (3.1)-(3.3) są ważne, gdy spełniony jest odpowiedni warunek linearyzacji (zgodnie z tabelką 5,1 oraz p. 2.1). W szczególności, warunek ten zapisuje się prosto:
(3.4) vKm ^ 0.1
dla sumowania „wzw” lub „pdw”.
3.1. Ocena wielkości zaburzeń
{t% }.Algorytm Romberga stosujemy za- zwyczaj dla niewielkich m, powiedzmy dla m < 15 (tzn. dla co najwyżej 215
~ 32760 węzłów). Wówczas z (3.3) otrzymujemy (16400 („wzw”), m < | 53 („pdw” i inne).
A więc sumowanie bardziej subtelnymi algorytmami gwarantuje małość zaburzeń (Ir^jJ ^ v • 53).
W przypadku sumowania „wzw” wielkość oszacowania dla K m może budzić pewien niepokój. Następujące dwie uwagi wskazują na często wystę- pujące sytuacje, kiedy zaburzenia mogą być akceptowane, mimo że charakterystyka K m nie jest mała.
U w a g i
(i) W przypadku symetrii zaokrągleń suma większej liczby pojedynczych błędów często zachowuje się jak suma niezależnych zmiennych losowych o wartości oczekiwanej zero i odchyleniu standardowym rzędu v. Dla więk- szych i oraz j jest więc bardzo prawdopodobne, że będzie spełniona nierów- ność
(lub że | t $J tylko nieznacznie będzie większe niż prawa strona tej nierównoś- ci).
Zatem dla j < i ^ m ^ 15 możemy oczekiwać (przy symetrii zaokrągleń), że będzie
1*211 < 130 („wzw”),
8 („pdw” i inne algorytmy).
(ii) Małość zaburzeń jryjj możemy rozumieć również w innym kontekś- cie, przez porównanie z błędami {<5/k(,)) obliczonych wartości {/k(,)], por.
p. 2.2. Równość
Tf' TU)
ljb — a
~2f~ l u w r
k= 0
wykazuje, że jeśli tylko zachodzi |/k(,) t JJJ ^ |ń/k(0|, to zaburzenia i f k tylko
nieznacznie powiększają niedokładność istniejącą w danych początkowych
Idealny algorytm obliczania wartości funkcji mógłby co najwyżej
Błędy zaokrągleń w algorytmie Romberga
37 gwarantować spełnienie |<5/k(,)| ^ v |/k(,)|. Aby sprostać takiej (najwyższej możli- wej) precyzji danych początkowych, należałoby zadbać o bardzo małą cha- rakterystykę kumulacji K m. Często jednak \óf£l)\ jest znacznie większe niż v\fk% co pozwala obniżyć wymagania „małości” zaburzeń {rfk}. □
3.2. Błąd względny obliczonej kwadratury. Oznaczmy przez 7}(,) dokładną wartość kwadratury dla obliczonych {/k(,)}:
(3.5) '
Z k= 0
oraz, przy założeniu 7}(i) / 0, wprowadźmy wytworzony błąd względny rów- nością
(3.6) = i(7}(£)— T)<i))/7}<i)|.
Z (2.35) i (2.36) otrzymamy równość przybliżoną
Zatem
S f =
r
' f a f p■+ŁQi+tT+xy] ,(i>-
(3.7) S f
i V[ r f ( r f - 2
j- 3)+(2/ + 3)], gdzie
(3-8) ry>= i"/3Mi/t0li/|I"^/*,i.
k = 0 k= 0
Mnożnik r f jest ^ 1, a równy 1 dokładnie wtedy, gdy wartości / k(i) są stałego znaku dla k = 0, 1,..., 2‘. Ostatecznie zapiszemy oszacowanie (3.91 q (0 < H 7? (2' - 1 + 4.5) + 2 ;> 3 ] („wzw”),
j "" (v [r f (i -f 4.5) -I- 2j + 3] („pdw” i inne).
3.3. Wykorzystanie regularności funkcji podcałkowej dla poprawienia osza- cowań. Dotychczas rozważaliśmy oszacowanie błędów dla zupełnie dowol- nych ciągów {/k(,)}k 2L0 danych początkowych. Tymczasem można oczekiwać, że są to wartości (na ciągu równoległych węzłów) funkcji / o pewnej klasie regularności.
A. Wykorzystamy najpierw założenie regularności funkcji / do „rozłado- wania” składników ekstrapolacyjnych { 0fk} zaburzeń na wszystkie dane Początkowe. Zgodnie z (2.39), w dotychczasowej analizie składniki te obcią- żały jedynie dane o wskaźnikach parzystych: {/ii1}-
Załóżmy mianowicie, że i jest dostatecznie duże, aby mogły zachodzić
równości przybliżone, por. ( 1.11):
38 J. Chojnacki, A. K iełb asi ńsk i
>* b
(3.10) 2‘~ 1
k = 0 2i~ 1Z Z I/2I - 1I = Pj,iy
/= 1 /= 1
(druga równość wynika ze wzoru prostokątów).
Wówczas dla wielkości
(3.11) A:= '% u f f l e v 2l
k = 0 1 = 0
z (3.10) wynika oszacowanie 2*~ 1
(3.12) M K v5 / l "/Sj. 2 ,1/ii’l = v 5 ,( i - ^ ) g . Przyjmując
(3.13) 6% = — sign(/k(i)),
otrzymujemy równość (3.14)
fc= o oraz
(3.15) \0%\ ^ Ą , 0j 1 - fc.i
Uwzględniając (1.11) oraz tablicę 2.2, otrzymujemy tablicę 3.1.
11. Tablica oszacowań \dj)
i0 0.000
1 0.333
2
0.963
3
1.333
4
1.451
5
1.491
6 1.503
7 1.506
Porównując (3.11) i (3.14), możemy zastąpić zaburzenia zaburzeniami
„wyrównanymi” Ofy, o lepszym oszacowaniu globalnym. Przy założeniu
(3.10) możemy więc skorygować oszacowania (2.37), (3.2) i dalsze, zmniejsza-
jąc je o 4.
Błędy zaokrągleń w algorytmie Romberga
39 B. Z regularności funkcji podcałkowej wynikają niekiedy zależności, po- zwalające na poprawienie oszacowań (3.9) błędu względnego obliczonej kwa- dratury. Dla i dostatecznie dużego zachodzą bowiem równości przybliżone
(3.16) Si? = C: = -r^— b - a i f/(x)rfz, Z (i):= S N (i)/2i~ ł ^ C (druga równość wynika ze wzoru prostokątów).
Jeśli liczba C jest wyraźnie różna od zera, to dla i dostatecznie dużego oba składniki sumy (2.11) są tego samego znaku, a ze wzrostem i zbliżają się do jC. Możemy to wykorzystać do innej niż dotychczas konstrukcji zaburzeń P(i)
b 0,k-Przedstawiając SN(i) w postaci, por. (5.2),
(3.17) SJV‘‘> = [ £ f i i l . U + ^ a + ■/,<?),
i = 1
(3.18) W 0! s: v p , ^ ’| 4 v^g>
oraz przyjmując założenie indukcyjne (3.i9) s ' r ‘> = ^7i T2£
k = 0
(3.20) przy (3.21)
1-4 ' 1}| ^ vź(l X), |jU(‘ X)| ^ vfll 1),
X<0) = ^ = F * = 0 , / i (1) = l ,
otrzymujemy z (2.11) (3.22)
gdzie (3.23)
Ą ° = i I " /" ’( ! + 4 °)(1 + V 0).
^ k = 0
/li0 = Ąl' l\ k = 2l, (3.24)
(3.25)
i j/f\ k = 2/ - 1, /i(i) == IJi + di,
nt = (sr V ,'~ 1 )+z(V(?)/(2Si?),
Z (,) = S N ® /? '1.
Prawdziwe są przy tym oszacowania
(3-26) l/L^I ^ vX(,), X(,) = max(X(,_1), ^ (,)),
(3-27) |/*(i)| £w = /7ł + 1,
(3-28)
40 J. Choj nacki, A. K iełb asi ńsk i
Przy założeniu
(3.29) S g - ‘»Z<'»>0, i> p , zachodzi
(3.30) n, t f ) ,
a przy założeniu
(3.31) S t 1] i ^ p ,
zachodzi
(3.32) 77, « ( / ? '- 11+ TO /2.
Zależności te pozwalają wyprowadzić zestaw oszacowań (przy założeniu (3.29) lub (3.30)), przedstawiony w tablicy 3.2.
3.2. Tablica wielkości XU),
i ^ pAlgorytm sumowania x<» p*0 (3.29) (3.30)
wzw 2‘~2 — 1 2‘~2+ 1 f 2' ~ 2 4- 2
gm, wzwfl2 max(l, p —2)
( i- p + 2)3
pdw
i- 2 22
Porównanie (3.22) z (2.12) daje równość przybliżoną
(3.33) + A
w której składnik może mieć (w przypadku niektórych algorytmów sumowania) istotnie lepsze oszacowanie niż składnik x{k z równości (2.13).
C. Reasumując, założenie regularności funkcji podcałkowej pozwala czę- sto na skorygowanie oszacowania zaburzeń (3.2) do postaci
(3.34) f<» = i 2‘" 1 + ^ + 3-5 (»wzw”>.
J ]i+ 2 / + 3.5 („pdw” i inne),
a oszacowania błędu względnego (3.9) do postaci (przy p = l w (2.29)) ( [ r • (2‘ - 2 - 0.5) + (2‘ - 2 + 2; + 4)] („wzw”),
(3.35) S f < v (T • (i — 0.5)+ (2j -f 4)] („pdw”), ( [ r • 1 -f (i + 2j + 2.5)] („gm, wzw fl2”), gdzie
(3.36) r:= ]\f(x)\dx/\]f(x)dx\.
a a
Błędy zaokrągleń w algorytmie Romberga
41 3.4. Użycie innych arytmetyk numerycznych. W rozdziale 2 rozważaliśmy zastosowanie arytmetyki binarną z najlepszą możliwą realizacją działań arytmetycznych. Nie wszystkie arytmetyki zmiennopozycyjne mają te włas- ności.
Obecnie założymy, że:
— realizacja działań całkowitych jest dokładna,
— realizacja mnożeń lub dzieleń liczb rzeczywistych (lub rzeczywistej i całkowitej) odbywa się zgodnie z regułą
fl(*Oy) = x O r ( i + e ) ,
O e
{*,/}, N < V.
(Dotyczy to również mnożenia (dzielenia) liczby rzeczywistej przez całkowitą potęgę liczby 2).
— realizacja dodawań i odejmowań liczb rzeczywistych odbywa się zgod- nie z regułą
(3.38) fl(*±.y) = x -(1+ s 1)±>'(1+£2)>
|£ ;|^ v , i = 1, 2.
Tę ostatnią relację nazwiemy złym sumowaniem.
Algorytmy gm i gm2, opisane w rozdz. 5, nie gwarantują w tym przypad- ku podanych w tablicy 5.1 oszacowań wielkości <p(,). Dla pozostałych algoryt- mów sumowania nadal jednak są prawdziwe oszacowania <p(l), a częściowo również ip }, {//$ z tablicy 5.1. Odpowiedniej modyfikacji muszą ulec zależnoś- ci (2.13) i (2.15) analizy błędu, a w konsekwencji również przedstawienie zaburzeń xftk w (2.35) oraz oszacowań. Zachodzi mianowicie dla obliczonej wartości 7}(,):
(3.39) Tf> = b — a Z I + t JU
k= 0 iTftl
(3.40)
T<‘>t m [zi0+ % ] + K}0+xł°], (3.41) fW = Lj [2‘- 1 + 6.3] + [3./ + 3]
[ i+ 6.3]+ [3 /+ 3 ]
(„wzw”),
(„pdw, wzwfl2”),
przy czym 9fk oznacza tu „wyrównany” składnik ekstrapolacyjny (przy założeniu regularności funkcji, por. p. 3.3):
(3-42) f t | ^ v 6 . 3 .
Porównanie z (3.34) wykazuje, że użycie gorszej arytmetyki pogarsza ogólne oszacowanie zaburzeń tylko nieznacznie.
4. Inne warianty algorytmu Romberga
4.0. Najprostszą modyfikacją algorytmu 1, często zalecaną (por. [5]),
42 J. Chojn acki, A. K iełb as iń sk i
Algorytm 2
Algorytm 2 różni się od algorytmu 1 jedynie zastąpieniem instrukcji (1.17) instrukcją
(4.1) S f: = Sf. 1 +(Sfl j 1)
W przypadku zastosowania opisanej w p. 2.0 arytmetyki można wykazać, że zachodzi
(4.2) T.b) ij b — a
k= o
(4.3) I t J'11 ^ v[<p(0+ i + 5],
przy czym udział „wyrównanego” składnika ekstrapolacyjnego §fk w tym oszacowaniu wynosi v-2. Zastosowanie sumowania pdw, gm, gm2, wzwfl2 gwarantuje więc oszacowanie kumulacji wielkością
(4.4) Tj* — 1+7 + 4 ^ K m — 2m + 4 .
W tym sensie algorytm 2 jest nieco lepszy od algorytmu 1 (dla tych algorytmów sumowania). Przy sumowaniu wzw przewaga algorytmu 2 staje się zupełnie nieistotna.
4.1. Powstaje pytanie, czy istnieje wariant algorytmu Romberga, dla którego charakterystyka fj0 kumulacji błędu zaokrągleń byłaby niezależna od wymiaru m zadania, j ^ i ^ w?
W przypadku arytmetyki opisanej w p. 2.0 i przy warunku linearyzacji typu
(4.5) v 4 ' < l ,
może to zapewnić następujący algorytm, wykorzystujący wariant (4.1) wzoru ekstrapolacyjnego, algorytm gm do obliczania SN(i), oraz poprawki typu G illa-M ollera do wszystkich pozostałych sumowań. Poprawki te tworzą tabli- cę trójkątną {ef}i = 0....mj= o , p o w s t a j ą c ą równocześnie z tablicą {Sf}.
Algorytm 3. Algorytm oblicza m +1 pierwszych wierszy tablic trójkątnych {Sf}{ef}i =
0 , . . . , m j = o , . . . , itakich, że
(4.6) T f):= (b -d )(S f + Ą )).
Algorytm 3
(4.7) S<00): = ( / r + / n / 2;
for i: = 1 to m do begin
2ł~ 1
(4.8) SN {i):= Z fzi— i > (gm) /= 1
(4.9) S f: = 1'/2 + SJV<i,/2‘;
Błędy zaokrągleń w algorytmie Romberga
43 (4.10) e®: = e% 1)/2 + (SN(0/2* i-(S 8 )-- S r ‘V2));
for j: = 1 to i do begin
(4.11) óf: - i) ;
S f: = S% 1 + d f;
p(i) • end end;
Można wykazać, por. [4], że zachodzi w tym przypadku równość
(4.12) + W
1 k=0
z oszacowaniem zaburzeń (przy „wyrównanym” składniku ekstrapolacyjnym)
(4.13) \x%\ ó - 8
oraz oszacowaniem błędu względnego, por. (3.35),
(4.14) ^ v [ / y 5 + 3].
U w a g a. Zastąpienie w instrukcji (4.8) algorytmu gm algorytmem gm2 pozwala złagodzić warunek (4.5) do
(4.15) v • 2* < 1 ,
przy równoczesnym wzroście oszacowań
(4.16) | t & |^ v 10,
(4.17) & P ś v t r - 7 + 3]. □
5. Algorytmy sumowania. Przypominamy tu ważniejsze fakty, dotyczące kilku algorytmów sumowania oraz ich implementacji.
5.0. Rozważane zadanie to obliczenie w arytmetyce fl (por. p. 2.0) sumy (5.1) s* = £ ah n = 2i~ \ / > 2,
/= i
dla danych liczb tej arytmetyki.
Każdy z przedstawionych w p. 5.1 algorytmów nadaje zmiennej rzeczywis- tej s, jako jej końcową wartość — obliczoną sumę. Analiza błędów zaokrąg- leń (por. [2], [3], [4], [6]) wykazuje, że istnieją zaburzenia {ę>J, {ipt}, takie że zachodzi równość
(5.2) s = Ż tf|(l + <Pi)= ż M l + iAiMl + ^o),
1= 1 /= 1
I
44 J. Chojn acki, A. K iełb asi ńsk i
oraz oszacowania
(S. \\(pi\ ś v<p(i), |^,| ^ vtf«\ / = 1,..., n, 1 |«Aol ^ vl?(o (<P (0 4= «A(i) + «Ao)•
W tablicy 5.1 zestawiamy charakterystyki kosztu i dokładności rozważanych algorytmów.
Charakterystyki dokładności (w postaci mnożników ę (i), f (l), ij/^) otrzymy- wane są przez linearyzację odpowiedniej funkcji błędu, są więc ważne jedynie wówczas, gdy spełniony jest odpowiedni warunek linearyzacji (por. p. 2.1), który zamieszczamy również w tabelce.
Tablica 5.1. Koszt i dokładność algorytmów sumowania (n = 2‘ *)
Algorytm Liczba sumowań
w fl
Liczba działań całkowitych
Obcią- żenie pamięci
Warunek linearyzacji
f. błędu
ę {i)tp> </'o)
wzw
n n2
vn < 0.1 n - 1(n/2)-l
n/2pdw
n Inlog2 n + 5 vlog2 n < 0.1 log2 n log2 n —l 1*
gm 4 n
n5
vn2 < 1 2*1* 1*
gm2
4 n n11
vn < 0.1 4.5*3.5* 1*
wzw fl2
n** n4 vn < 1
21 1
* — oznacza pozycje, które nie są prawdziwe dla arytmetyki ze złym zaokrągleniem sumy, p. (3.38).
** — oznacza sumowanie w rejestrze podwójnej precyzji.
Tablica zawiera oszacowania globalne. Dokładniejsza analiza wykazuje m.in., że
(5.4)
(p{2) = 1, iJ/{2) = 0
<p(3) = 2 , <p{4) = 3
i£(3) = l, i?(4) = 2
dla gm i gm 2, dla gm 2.
5.1. Algorytmy. Algorytmy A, B, C, D korzystają wyłącznie z arytmetyki fl oraz arytmetyki liczb całkowitych.
Zmienne typu real: s, z, u, v, x, y real array: ss[l:lo g n ] integer: p, k, l, q, j
Algorytm E korzysta ponadto z rejestru sumowania w podwójnej precyzji, oznaczanego tu jako zmienna
real double spp
Błędy zaokrągleń w algorytmie Romberga
45 A. Algorytm wzw („wyraz za wyrazem”)
s: — 0;
for /:= 1 to n do si—s + ap,
B. Algorytm pdw („po dwa”)
P - — 0 ;