Błędy zaokrągleń w algorytmie Romberga(Praca wpłynęła do Redakcji 1985.07.04)

Pełen tekst

(1)

ROCZNIKI POLSKIEGO TOWARZYSTWA MATEMATYCZNEGO Seria III: MATEMATYKA STOSOWANA XXIX (1987)

J

a n u s z

C

h o j n a c k i

, A

n d r z e j

K

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-

(2)

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

+

X

E,

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).

(3)

Błędy zaokrągleń w algorytmie Romberga

27 Kwadratury te tworzą tablicę trójkątną

t u )

73'01 T

q

1) 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 @

q

a 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.

(4)

28 J. Chojn acki, A. K iełb as iń sk i

gdzie

( 1 . 10 )

Po,k ~

1>

P j - i,k/2

Pj-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)).

(5)

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ą.

a

Ponieważ 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

(6)

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.

(7)

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)( \ +

q

0), 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.

(8)

32 J. Chojn acki, A. Kieł basi ńsk i

Przyjmując założenie indukcyjne

2

i_

1

S (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) + +

(9)

Błędy zaokrągleń w algorytmie Romberga

33

(2.19) 1,2/'

e(« - D l

k = 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

(10)

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

(11)

Błędy zaokrągleń w algorytmie Romberga

35 T

w ie r d z e n ie

1. (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,l

j

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).

(12)

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)

lj

b — 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

(13)

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):

(14)

38 J. Chojnacki, A. K iełb asi ńsk i

>* b

(3.10) 2‘~ 1

k = 0 2i~ 1

Z 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)

i

0 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.

(15)

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)

(16)

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 ^ p

Algorytm 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 2

2

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

(17)

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]),

(18)

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 , . . . , i

takich, ż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‘;

(19)

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

(20)

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 n

2

vn < 0.1 n - 1

(n/2)-l

n/2

pdw

n In

log2 n + 5 vlog2 n < 0.1 log2 n log2 n —l 1*

gm 4 n

n

5

vn2 < 1 2*

1* 1*

gm2

4 n n

11

vn < 0.1 4.5*

3.5* 1*

wzw fl2

n** n

4 vn < 1

2

1 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

(21)

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 ;

for /:= 1 to n do begin

p: = P + l ; q'-= / - l ; s:= at;

•RET: k : = q ^ 2;

if 2*k = q then ss [p ]: = s else begin

q:=k\ p: = p - 1; s : = s s [ p ] + s;

go to RET end; end

C. Algorytm gm (Gilla-Mollera, por. [4]) if n = 2 then s: = + a2

else begin

s:= z: = 0;

for /:= 1 to n do begin

u: = at; u: = s + w;

z: = z + (w — (u — s)); s:=i>

end; s : = s + z end;

D. Algorytm gm2 (dwukrotny Gilla-Mollera, por. [3]) if n = 2 then s: = flj a2

else

if n = 4 then s:= (ax+ a 2) + (a3 + n4) else begin

p:=entier (sqrt(n));

if p*p ^ n then p := entier(sqrt (n 2));

q:= n-j- p; s:= z: — 0 ;

(22)

46 J. Chojn acki, A. Kieł basi ńsk i

for k:= 1 to ą do begin

j: = k*p\

if p = 2 then begin x:= aJ+l\ y:= aj+2 end else

begin

for /: = (k — l)*p + 1 to j do begin

u:= at; v\— x + u;

y:= y + (u -(v-x));

s:~ v end

end; u: = x + y ; v:= s + w;

z: = z + (u—(v — s));

s:= v end; s: = s + z

E. Algorytm wzw fl2 (por. [2], [6]) spp:= 0;

for /:= 1 to n do [spp:= spp + a j -fl2;

s:=spp;

Prace cytowane

[1] J. C hojnacki, Metoda Romberga, praca magisterska, Stud. Zaoczne Wydz. Mat., Inf., Mech. UW 1982.

[2] J. i M. Jankowscy, Przegląd metod i algorytmów numerycznych, Warszawa, WNT 1981.

[3] M. Jankowski, A. Smoktunowicz, H. Woźniakowski, A notę on floating-point

summation of very many terms, Elektron. Informationsverarb. u. Kybernetik 19 (1983), 9,

435-440.

[4] A. Kiełbasiński, Algorytm sumowania z poprawkami i niektóre jego zastosowania, Mat.

Stos. I (1973), 22-41.

[5] J. Stoer, Wstęp do metod numerycznych, Warszawa, PWN 1979.

[6] J. H. Wilkinson, Błędy zaokrągleń w procesach algebraicznych, Warszawa, PWN 1967.

Obraz

Updating...

Cytaty

Updating...

Powiązane tematy :