• Nie Znaleziono Wyników

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

N/A
N/A
Protected

Academic year: 2021

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

Copied!
22
0
0

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.

Cytaty

Powiązane dokumenty

Nierówność (2.1) jest uogólnieniem nierówności Bhattacharya rzędu k. Na przykład, gęstość rozkładu jednostajnego.. gdy chcemy podać dolne ograniczenie typu RC dla funkcji

Rozpatrywany model pojawia się wtedy, gdy szum urządzenia pomiarowego zniekształca stan układu w chwili pomiaru i wpływa na stan układu w chwili następnej.. Ponadto podajemy

Niech H będzie dwuwymiarową dystrybuantą rozkładu dodatnio kwa- drantowo zależnego o brzegowych rozkładach jednostajnych na odcinku [0,1]. Okazuje się, że jest

Celem pracy jest wyjaśnienie, w jakim sensie można mówić o rozwiązywaniu zadania źle postawionego i pokazanie, że stosując odpowiednią metodę udaje się zadanie

1» Praca R.F.Greena dotyczy uproszczonego przypadku, w którym osobniki (zwierzęta , rośliny) nie konkurują ze sobą o zasoby środowiska (np. Przyjęcie takiego uproszczenia

Nie ma sensu mówić o zbieżności występujących powyżej szeregów zmiennej h, ponieważ nie ma wzoru Taylora w przestrzeni D'(Rn) (por. Pokażemy istnienie

Najwybitniejszym naukowcem w początkowym okresie historii tej dziedzi- ny był niewątpliwie Vilfredo Federico Damaso Pareto (1848-1923). są uważane za podstawę

Gleser, On the asymptotic theory of fixed size sequential confidence bounds for linear regression parameters, Ann. Gołdys, Stałoprecyzyjna estymacja średniej