WrocławNiektóre zagadnienia związane z absolutną stabilnością jawnych metod jednokrokowych

M i e c z y s ł a w S z y s z k o w ic z

Wrocław

Niektóre zagadnienia związane z absolutną stabilnością jawnych metod jednokrokowych

(Praca wpinęła do R edakcji 4.10.1983)

1. WSTĘP

W pracy tej rozważa się Jawne metody Jednokrokowe rozwiązywania za- gadnienia początkowego

(1.1) y' = f(x,y), y(xQ) dane.

Zagadnienie (1.1) może być Jednym równaniem różniczkowym bądź ukła- dem równań różniczkowych. Zakłada się, że rozwiązanie powyższego zagadnienia istnieje i Jest Jednoznaczne. Rozwiązanie to Jest dalej oznaczane przez y(x).

Metoda Jednokrokowa Jest najczęściej określona przez swoją funk- cję przyrostu $f(x,y,h), gdzie przy danych xn, yn, h^ przybliżone rozwiązanie zagadnienia (1.1) w punkcie *n+hn otrzymuje się na pod- stawie wzoru

[163]

0.2) yn+1 = y„+hnffcxn,yn,hn).

Wartość Yq = y(x0) Jest zadana w warunku początkowym w (1.1). Wiel- kość hn nazywa się krokiem całkowania metody jednokrokowej i naj- częściej hn zmienia się w trakcie prowadzenia obliczeń ([²]). W dal- szej części tej pracy funkcja (x,y,h) jest oznaczana przez $, a metoda jednokrokowa jest utożsamiana z funkcją Opuszcza się rów- nież wskaźnik n przy x oraz h. Tak więc, przy przyjętych oznacze- niach, wzór (1.2) ma postać

0-3) yn+i - yn+h?.

W pracy niniejszej o metodzie jednokrokowej r zakłada się, że zastosowana do zagadnienia początkowego

(1.4) y ' = A y, y(x0) dane, AeC,

daje następujący związek między wartościami rozwiązania przybliżo- nego yn i yn+1

01.5) yn+1 = w(z)yn ,

gdzie z = hr\ oraz

2 _ ni .

(1.6) w(z) = 1 + z + —— +■ ... + —— + y d. — . 2 * T)l p - i=p+1 1 1 i I11

Liczba m jest nazywana stopniem metody jednokrokowej. Współczynni- ki d^ w (1.6) zależą od parametrów metody jednokrokowej<§. W lite- raturze (np. [³]) wielomian w(z) jest nazywany wielomianem stabil- ności.

DEFINICJA 1.1. Wielomian stabilności w(z) ma rząd zgodności q, jeżeli spełnia następujące związki:

d^t-w(z) = 1 przy z-»0, j = 0, 1, ..., q. ri dzJ

Zachodzi następujące twierdzenie ([3]).

TWIERDZENIE 1.2. Metodzie jednokrokowej $ rzędu q odpowiada/

wielomian stabilności w(z) , którego rząd zgodności wynosi q.j—j

Przy stosowaniu metody $ do rozwiązywania zagadnienia (1.1) jest wymagane, aby błąd globalny metody, tj. £n = yn-y(xj » przy ustalonej długości kroku całkowania był ograniczony, gdy n —^*00 .

Badanie tej własności w przypadku dowolnego zadania (1.1) jest na ogół trudne. Z tego też powodu używa się zwykle zadania testowego

(1.4). Ma to następujące uzasadnienie.

Rozwiązanie zagadnienia (1.1) z funkcją f(x,y) różniczkowaIną względem y, lokalnie zachowuje się jak rozwiązanie zagadnienia li- niowego

y' * Ay,

gdzie macierz A Jest Jakobianem ^ • Jeżeli macierz A ma s (s ok-Sf reśla liczbę równań w (1.1)) różnych wartości własnych -*2, •••*

Ag, to istnieje taka macierz nieosobliwa P, że funkcja z « P y -1 spełnia s niezależnych równań

z' = Dz,

gdzie D = diag(A1tA2, ..., Ag). Tak więc, ograniczoność błędu glo- balnego yn-y(xn) dla n-*oo Jest równoważna ograniczoności błędu globalnego zn-z(xn) . Liczba AeC w (1.4) reprezentuje dowolną war- tość własną --- •df

ay

DEFINICJA 1.2. Metoda $ zastosowana do zadania (1.4) Jest ab- solutnie stabilna dla danej liczby hAeC, gdy ciąg otrzymany w (1.3) yn-»0 przy n — ► oo •

Z daną metodą Jednokrokową związany Jest obszar absolutnej sta- bilności, który można zdefiniować następująco.

DEFINICJA 1.3. Podzbiór KcC nazywa się obszarem absolutnej sta- bilności metody $ , Jeżeli metoda $ Jest absolutnie stabilna dla każdej wartości hAeK.

\

DEFINICJA 1.4. Przecięcie obszaru absolutnej stabilności K z osią rzeczywistą nazywa się przedziałem absolutnej stabilności me- tody. Przedział ten Jest w tej pracy oznaczany przez L CL » KnR).

Powyższe definicje można znaleźć np. w pracach [3] , [7] , fi 5*] • Szczególnie dużo uwagi metodom Jednokrokowym i ich własnościom Jest poświęcone w [3] •

W pracy niniejszej są podane przykłady metod jednokrokowych.

Przedstawione są metody jednokrokowe Bobkowa o rozszerzonych obsza- rach absolutnej stabilności. Zbadane jest zachowanie się przedzia- łów absolutnej stabilności metod jednokrokowych po ekstrapolacji Richardsona. Zaproponowany jest taki dobór parametrów metody $ , aby w przypadku równania (1.4), otrzymać po ekstrapolacji Richard-

sona najwyższy rząd dokładności. Przedstawiona jest metoda rozpoz- nawania charakteru zależności kroku całkowania, tzn. rozstrzyganie, czy wielkość kroku całkowania zależy od żądanej dokładności, czy też od długości przedziału absolutnej stabilności metodyk. Podane są przykłady z obliczeń na maszynie cyfrowej.

2. PRZYKŁADY JAWNYCH METOD JEDNOKROKOWYCH

Jednymi z bardziej znanych metod jednokrokowych, są metody Rungego- -Kutty (RK)• Ogólnie m-stopniowa jawna metoda RK może być zapisana w postaci

m

yn+1 " V hi- $ - I °rkr>

r=1

gdzie

ar r-1

^ri = ****

i=1

Niech p*(m) oznacza najwyższy rząd metody RK m-stopniowej. Zacho- dzą następujące związki ([⁷]):

p*(m) = m, m « 1, 2, 3, 4,

pM(m) = m-1, m = 5, 6, 7,

(².¹)

pK(m) = m-2, m = 8,9,

PK(m) ^ m-2, m = 10,11, ...

W przypadku p-stopniowej metody RK rzędu p (p = 1, 2, 3, 4) otrzy- muje się ([7j):

Tabela 2.1

Rząd metody Wielomian w (z) Przedział L

p = 1 1+z (-2.0,0)

p = 2 1+Z+— z2 (-2.0,0)

2 !

p = 3 i+z+5_+z_ 2 _3 (-2.51,0) 2! 3!

2 3 4

p = 4 i+z+5_+5_+z_ (-2.78,0)

2! 3! 4!

Dla metod RK rzędu p z p y 4, na podstawie (2.1) oraz twier-

dzenia 1.2, otrzymuje się następujący wielomian stabilności w (z):

Tak więc, wielomian stabilności zależy od metody RK. W szczególnoś- ci, można tak dobrać parametry metody RK rzędu p, aby wielomian w(z) miał rząd zgodności m, tzn. aby był postaci

Taki wielomian stabilności uzyskuje się zawsze w przypadku jawnej metody Taylora. W metodzie Taylora funkcja $ jest zadana następu-

Metoda Taylora jest rzadko używana w obliczeniach na maszynie cyf- rowej ze względu na konieczność wyznaczania wyższych pochodnych

stabilności metody Taylora dla m = 1, 2, ..., 7.

Parametry metody RK można również dobrać tak, aby wielomian ok- reślony w (2.2) miał możliwie najdłuższy przedział, na którym

|w(z)| <1. Innymi słowy, aby przedział absolutnej stabilności L me- tody RK był najdłuższy.

(2,2) w(z) = 1+Z+...+ — +ZP p*

(2.3) w(z) = r — r~L i!

funkcji y(x). W pracy [a] przedstawione są obszary K absolutnej

PRZYKŁAD 2.1. Lawson [8] zaproponował następującą rodzinę metod RK rzędu 5 i 6-stopniowych:

21 1 2

41 ii 1

16

21 ⁴¹ - 166 — 166 326

"43 16 + 126 - if + 125 246 9 TC

1 4

7 * 1926

7 1 1926i - ? 3846'

7 12

~ 7 8 7

JZ90 0 .22

90 12

90 32

90 7

90

Dla metod z tej rodziny wielomian (2.2) jest postaci

w (z) = y + 366 — i=0 i! 61

gdzie 6 jest parametrem. Lawson [8] przyjął 0 * co daje prze- dział absolutnej stabilności L = (-5.62,0), Ze względu na kształt obszaru K zrezygnował on z tej wartości 6 , przy której przedział L jest najdłuższy.

Obok wymienionych już metod RK oraz metod Taylora, można podać jeszcze jeden przykład jawnych metod jednokrokowych - metody jedno- krokowe Bobkowa (RKB)• Opis tych metod można znaleźć w oryginalnej pracy Bobkowa [i], w książkach [5,6] lub też w pracy [i 3], gdzie podany jest algorytm otrzymywania metod RKB.

Tutaj podana jest jedynie idea konstrukcji metod RKB. Metodę RKB rzędu p otrzymuje się na podstawie żądania, aby rozwinięcie

<1

względem potęg h wyrażenia y(x+h)-y(x) oraz wyrażenia h A.y'(x+oc,h) i=0 1 1 było zgodne włącznie do wyrazów zawierających łr . Otrzymuje się następujący układ równań na A^ orazoc^:i

q q

(2• 4) A^ = 1, ^ Ai ^i = j+1 ( ^ = ^ ^» * * * * P*”^) •

i=0 i=0

Często {bc } nazywa się wagami, a (A.J - węzłami metody. Z te- orii numerycznego całkowania wynika, że jeżeli p-^2q+1, to układ powyższy ma rozwiązanie. Nieznane wartości y(x+0c^h) można znów ob- liczyć tak jak wartość y(x+h). Wystarczy teraz uzyskać je metodą rzędu p-1• Metody RKB określa się jako metody rekurencyjnego pod- wyższania rzędu.

Do zapisywania metod jednokrokowych Bobkowa wygodnie jest użyć następujących oznaczeń;

yn+oc - przybliżenie wartości y(xn+och),

yn+oę - przybliżenie wartości y(xn+och) z zaznaczeniem, że lokal-Ir ny błąd tego przybliżenia jest 0(h ) / Ic

i analogicznie

£n+ix * f(xn+Cch- W oraz ^n+oc * f(xn+«h* ^ +oc)‘

PRZYKŁAD 2.2. Metoda RKB drugiego rzędu i 2-stopniowa

y«+1 “ yn+h£n’

^ +1 - + |(fn+fn+l)’ □

Przy konstrukcji metod Jednokrokowych Bobkowa Jest duża dowol- ność stosowania wzorów dających rozwiązanie w węzłach {oc^} z żąda- ną dokładnością. Metody Bobkowa mają tę własność, że każdy wzór wchodzących w ich skład określa wartość rozwiązania przybliżonego w pewnym punkcie przedziału [x,x+h]. Takiej własności nie mają np.

metody RK, gdzie wartości liczone w punktach przedziału [x,x+h] nie są na ogół wartościami rozwiązań zagadnienia początkowego w tych punktach. Ta szczególna własność metod RKB daje się wykorzystać w metodach doboru kroku całkowania ([13])*

3. METODY JEDNOKROKOWE BOBKOWA 0 ZWIĘKSZONYCH PRZEDZIAŁACH ABSOLUTNEJ STABILNOŚCI

Długość przedziału absolutnej stabilności metody JednokrokoweJ $ Jest określona przez zachowanie się odpowiadającego jej wielomianu w(z) , tj. wielomianu określonego w (1.6). W przypadku metod RK rzę- du p = 1, 2, 3 i 4 oraz p-stopniowych, przedziały te są określone Jednoznacznie (tabela 2.1). Dla metod RK wyższych rzędów poprzez odpowiedni dobór parametrów metody RK można uzyskać zwiększenie długości przedziału absolutnej stabilności (przykład 2.1, [8,9]).

Doborem współczynników d^ (i = p+1, p+2, ••*, m) w wielomianie w(z) z (1.6), tak, aby |w(z)|^1 na możliwie najdłuższym przedzia- le, zajmował się Riha [lo]. Użył on do tego celu zmodyfikowanego algorytmu Remeza.

Van der Houwen w swojej książce [3] podał wartości współczynni- ków d^ oraz końce przedziałów absolutnej stabilności. Wyznaczył on

również takie wartości d^, przy których wielomian w(z) ma najdłuż- szy przedział na osi urojonej, na którym to przedziale |w(z;| ^ 1.

Poza tym Van der Houwen podał kilka sposobów otrzymywania wartości d^, przy których |w(z)| ^ 1 na najdłuższym przedziale.

Konstrukcja metod jednokrokowych o możliwie najdłuższym prze- dziale absolutnej stabilności jest podyktowana względami praktycz- nymi. Zdarza się, że rozwiązując zagadnienie początkowe (1.1), stwierdzi się następujący fakt.

Koszt rozwiązania zagadnienia, mierzony liczbą obliczeń war- tości funkcji f(x,y) (oznaczenie: [f]), nie zależy od żądanej dok- ładności £c[£q, Krok całkowania jest zdeterminowany wielkoś- cią przedziału absolutnej stabilności metody (;1l]). Obrazuje to ry- sunek 3.1.

Rys. 3.1. Zależność [f] od £ . Dla £ ^ £ >j krok całkowania jest zdeterminowany wielkością przedziału L

W takiej sytuacji zastosowanie metody 3? o długim przedziale absolutnej stabilności urealnia zależność kosztów od wymaganej do- kładności. W paragrafie 6 tej pracy podany jest sposób automatycz- nego wykrywania sytuacji przedstawionej na rysunku 3.1.

Tutaj przedstawione są metody jednokrokowe Bobkowa rzędu p =

= 1, 2, 3 oraz 4 i ro-stopniowe, gdzie m = p+1, co pozwala uzyskać metody o większym przedziale absolutnej stabilności niż w przypad-

ku metod RK rzędu p i p-stopniowych. Podane są przykłady takich metod, a w niektórych przypadkach całe ich klasy.

Absolutna stabilność metod pierwszego rzędu

Na podstawie wzoru (2.4) dla metody RKB pierwszego rzędu otrzymuje się równanie

q

Z Ai = 1 * i=0

' Parametry oc^ (i = 0, 1, ..., q) mogą być dowolne. Chcąc Jednak ko- rzystać tylko z wartości na przedziale [x,x+h], wprowadza się ogra- niczenie

(3*1) o ^ oci < 1 (i = o, 1, ..., q) •

Do wyznaczenia metody RKB pierwszego rzędu 2-stopniowej można przy- jąć dowolną wartość q. Otrzymuje się wtedy następujący wzór:

q

yn+i * v h E Aifn+Cc. •

i=0 1

gdzie wartość yn+cc oblicza się ze wzorów

yn+C(l - yn+orihfn Ci = 0, 1...q).

Tak określona metoda jest 2-stopniowa i między wartościami yn i yn+1 dla zagadnienia (1.4) zachodzi związek

/ q q

yn+i - (1łhA I! Ai+ (hA> E Aiai)v

q q

Przyjmując a = ^ ^^i oraz korzystając z tego, że ^ A^ = 1,

i=0 i=0

można określić wielomian zmiennej z

w (z) = az +z+1p

Znany jest następujący fakt ([lo]):

W przypadku gdy rząd metody p = 1, wielomian w(z), dla którego

|w(z)|<$ 1 na najdłuższym przedziale, otrzymuje się poprzez prze- kształcenie nHtego wielomianu Czebyszewa pierwszego rodzaju Tm w następujący sposób:

Najdłuższy przedział, na którym |w(z)| ^ 1, wynosi wtedy [-2m2,o].j-j Na mocy powyższego, wielomian

wQ(z) = az +z+1,p

przy a = -g- , ma najdłuższy przedział [-8.0,o] , na którym |w (z)|^ 1

Do realizacji metody przyjmuje się takie wartości A^, oc^, aby było |w (z)| < 1, co Jest spełnione, gdy a > g- .i

Przyjęcie q = 0 i ustalenie wartości a (oCq = 0) w sposób Jed- noznaczny określa metodę Jednokrokową

^ti+ocq “ yn+hoc0fn*

(3.2)

yn+i - yn+hW 0-

Jest to metoda, która wymaga dwóch obliczeń wartości funkcji f. ¥ przypadku q = 1, ocq = 0, dostaje się metodę

W , = yn+hV n >

(3.3)

yn+1 " yn+h (A0 f n+A1 f n+oc^ *

Metoda ta również wymaga dwóch obliczeń wartości funkcji f. Ustale- nie wartości a w tym przypadku pozwala otrzymać całą klasę metod, dla których musi być spełnione Aq+A^ = 1 oraz A^oc>| = a* Dla innych wartości q, wzrasta liczba obliczeń wartości funkcji f. Tak więc metoda określona w (3.3) Jest najbardziej interesująca z możliwych do otrzymania metod 1 rzędu 2-stopniowych.

Na rysunkach 3.2 i 3.3 przedstawiono obszary absolutnej stabil"

ności metody RKB pierwszego rzędu 2-stopniowej. Obszary są symetry- czne względem osi 0X.

1 1 Rys. 3«2. Obszary absolutnej stabilności przy a*^-, a - y •

Absolutna stabilność metod drugiego rzędu

Parametry metod drugiego rzędu muszą spełniać następujące równania:

q q

Z Ai * 1- E V i ■

i-0 i-0

oraz przyjęte ograniczenie (3.1).

Przykładem metody Bobkowa drugiego rzędu 3-stopniowej Jest me- toda określona wzorami

w - yn + f hfn«

(3-4J W - V * hfn+*/2'

yn+1 - V h#Ofnł U 1fn«'

gdzie przyjęto oc = o<^ oraz A^oc ® •A

Dla metody (3.4) otrzymuje się wielomian w (z) postaci

-3 2

w (z ) = ^OC — + — + ^Z + 1 . 4 2

Na podstawie badania przebiegu wartości wielomianu w^c’z) otrzymuje się, źe przy oc * -^ęA wielomian wK(z) spełnia warunek I w^(z) | 4 1 na najdłuższym przedziale i przedział ten wynosi [-6.26,0]. Parametr a jednoznacznie określa metodę (3.4). Do realizacji metody (3.4) można przyjmować oc z przedziału Dla oc = ^ obszar absolutnej stabilności metody przedstawiono na rysunku 3.4.

<4

Rys. 3.4. Obszar absolutnej stabilności metody (3.4) z oc= • Obszar jest symetryczny względem osi 0X

Absolutna stabilność metod trzeciego rzędu

Bobkow [5] podał następującą metodę trzeciego rzędu 4-stopniową:

yn+1/4 = yn + T fn*

yn+1/2 ” yn + T fn+1/4’

(3.5;

yn4-1 " V hfn+1/2>

yn+1 " ynł|(fn+4fn+1/2+fn+l)-

Metodzie tej odpowiada następujący wielomian w(z):

Na podstawie (3.5) można utworzyć rodzinę metod trzeciego rzędu.

W tym celu wystarczy wprowadzić parametr a i otrzymuje się:

yn+a/4 = yn + t hfn*

yn+a/2 = yn + T hfn+a/4*

(3.6)

yn+1 * yn+ahfn+a/2-

yn+1 - ynłh(A0fn+A1fn+a/2+A2fn+a):

Po zastosowaniu tych wzorów do zadania testowego (1.4) otrzymuje się związek

yn+1 = 0 + ( V A1+A2) h U (A1 t + A2a) (hA)2+

2 2 3

+ (A1 “8 + A2 ^2") (hA)5+A2 ^8" (h*/0yn*

(

Z warunków na rząd metody (3.6) dostaje się następujące równania:

Aq+A^+A2 1=1 1 i

A1 2 +A 2 a “ 2 »

A a? +A ai _ JL 1 8 2 2 ^“ 6 ^*

Dodatkowo wprowadza się parametr b i do tego układu równań dołącza się jeszcze jedno równanie

2 8 =

Rozwiązania powyższego układu równań można wyrazić w zależności od parametru b;

^ ■ 1- W Ai - (1 - ^ r ) i- . =

■i

oraz przy b ^ ^ otrzymuje się dwa rozwiązania na a

Wielomian zmiennej z = hA związany z metodą (3.6) ma postać

4 z^ z²

w^ (z) = bz + —g- + + 2+1 •

Poprzez ustalenie wartości parametru b określa się wartość a (przyj- muje się a.) lub a2) i w ten sposób otrzymuje się konkretną metodę postaci (3.6).

Na podstawie badań przebiegu wielomianu wb(z) otrzymuje się, że przy b^ "55 ([i'll) wielomian ma najdłuższy przedział, na którym

|wb(z)| ( 1 i przedział ten wynosi [-6.03,0].

Absolutna stabilność metod czwartego rzędu r ***

Bobkow w [6J przedstawił metodę czwartego rzędu 5-stopniową nastę- pującej postaci:

yn+1/3 " yn + ł f n«

yn+2/3 - yn + 3 h fn+1/3>

C5-8) yn+1 - yn + ł ( V 3 fn+2/3) >

yn+1/2 " yn + |(2V 3fn+2/3-fn+l) >

yn+1 ' yn + lK f n+4fn+1/2+fn+l)*

Wielomian w (z) dla metody (3.8) ma postać

5 4 3 2

, N z z z z j,

w (z) = - -72 + 24 + 5 + ~jT + z+^ •

W pracy [5] Bobkow podał klasę metod jednokrokowych rzędu 4 i 5*

-stopniowych

yn+l/l yn + 1 ^n*

„ . h (6-1 * , 1_ *

yn+1/3 yn 3 1 6 n 6 n+1/l) *

(3*9) yn+1/2 = yn + (fn+3fn+1/3)’

yn+1 = yn + T (fn"3fn+1/3+4fn+1/2)’

yn+1 “ yn + T (fn+4fn+1/2+fn+1^>

gdzie 1 jest parametrem (1? 0) i wartość 1 nie wpływa na postać wielomianu w(z):

Poniżej przedstawione są obszary absolutnej stabilności metody (3.8) i (3.9). Na rysunku 3.5 widać, że metoda (3.9) ma znacznie większy obszar absolutnej stabilności niż metoda (3.8), a zatem z praktycz- nego punktu widzenia lepiej jest używać metody (3.9). Najdłuższy przedział absolutnej stabilności dla metody rzędu 4 5-stopniowej wynosi C-6.06,0) ([11 ]).

Rys. 3.5. Obszary absolutnej stabilności metod (3.8) i (3.9)

Wyznaczanie obszarów absolutnej stabilności

Zadanie wyznaczenia obszaru absolutnej stabilności metody jednokro- kowej <£ polega na znalezieniu takiego obszaru KcC, aby

gdzie w(z) jest wielomianem określonym w (1.6).

W celu wyznaczenia obszaru K rozwiązuje się następujące równa- nie:

u

Z € K

w(z) = e ,it

gdzie t Jest ustaloną wartością rzeczywistą* Innymi słowy szuka się przeciwobrazu okręgu o promieniu 1. Do rozwiązania powyższego rów^- nania stosuje się metodę Newtona wyznaczania zera funkcji. Oznacza- jąc FCz,t) ss w(z)-e^ oraz przyjmując liczbę N i t^ » k (k ** 1, ' 2, ..., N), gdzie m Jest stopniem wielomianu w(z), a liczba N ok- reśla ilość punktów brzegu obszaru K, wykonuje się wskazane opera- cje związane z realizacją metody Newtona:

Ck = 1, 2, ..., N) (1 = 0, 1, ...)

Ponieważ brzeg obszaru K przechodzi przez początek układu współrzęd- nych, więc Jako przybliżenie z® bierze się punkt 0. Zbiór punktów {o, z^, ..*, Zjj} wyznacza brzeg K.

4. ABSOLUTNA STABILNOŚĆ METOD RUNGEGO-KUTTY PO ZASTOSOWANIU EKSTRAPOLACJI RICHARDSONA

W tej części pracy są rozważane metody RK rzędu p, gdzie p = 1, 2, 3 i 4. Zakłada się, że wyznaczone są daną metodą RK dwa rozwiązania yn+/j oraz Yn+-j» dla których na podstawie (1.3) i (1.5) spełniona Jest równość

.1+1 F^zk » V F *^zk>tk^

yn+i ■ A z W

Zwykle mając wartości yn+1 oraz yn+1 stosuje się ekstrapolację Ri- chardsona i wyznacza się wartość

(4 1) - y +

^ ' ^n+1 “ yn+1 2p ^ *

gdzie p jest rzędem metody. Tak więc stosuje się metodę Rungego- -Kutty-Richardsona (RKR), która daje rozwiązanie y^+1 i rząd meto- dy RKR jest p+1. Rozwiązanie £ * 1 może być traktowane jako rozwią- zanie otrzymane nową metodą jednokrokową RKR, której własności ab- solutnej stabilności są inne od metody RK. Na podstawie (4.1) otrzy- muje się

oraz związek

y£+1 - w*(z)yn.

W oznaczeniach przyjęto dla i 4 n*

Dla metody RKR przedział absolutnej stabilności jest określony jako zbiór tych z, że |w*(z)| < 1.

Niech Yn+>| oznacza teraz rozwiązanie otrzymane z krokiem h/k poprzez k-krotne zastosowanie metody RK, gdzie k = 2, 3, ..., 6.

Takie podziały kroku całkowania można stosować w praktyce w celu doboru kroku. Na długich przedziałach całkowania podziały z k = 3, 4, 5 są nieraz mniej kosztowne od podziału przez 2 ([l3]). Koszt mie- rzy się wielkością [f ] •

Analogicznie do (4.1) otrzymuje się

^ . yn+1~yn+1 ^ ^ r

^+1 “ yn+1 kp ^ * 53 2» 3» •••» 6,

oraz wielomian w*(z)

(4.3) w*(z) = w^f-J] + V(— , k = 2, 3, ..., 6.

VK' kp~1

Wyznaczenie wielomianów w*(z) oraz przedziałów, na których

< 1, jest w zasadzie wykonalne tylko za pomocą maszyny cyf- rowej. W celu wyznaczenia wielomianów określonych w (4.3) wykorzys- tuje się specjalny system automatycznych obliczeń analitycznych na wielomianach [14] . Następnie otrzymane wielomiany tablicuje się, aby móc wyznaczyć przedział absolutnej stabilności metody RKR.

Dla metody Eulera (tj. metody RK pierwszego rzędu) dostaje się następujące wielomiany w*(z) dla k = 2, 3, ..., 6:

w* (z) = + z+1,

= f § + 4 + z+1 -

= 192 + 12 + 2 + 2+1 *

w*(z) = z 5 4 z z_ z 3 2 ^ 2500 + 100 + 10 + 2 + 2+1 *

6 5 4 3 2

*, . z . z-' , z , z-'' z X, w (z) - 5qqqo + 1080 72 9 2

Dla metod wyższych rzędów wielomiany są bardziej złożone i najwyż*

szy stopień wielomianu w^cz) wynosi 24.

Przedziały absolutnej stabilności metod RKR są przedstawione poniżej w tabeli 4.1 w postaci (a,0), gdzie podane są tylko war- tości a i pominięto znak minus.

Tabela 4.1

p k 2 3 4 5 6

1 2.0 3.0 4.0 5.0 6.0

2 5.15 7.15 9.15 11.15 ^13.10

3 4.10 6.50 8.10 9.90 11.70

4 6.45 9.15 10.20 12.10 14.40

Z tabeli 4.1 można wysnuć następujące wnioski:

1° dla p = 1, 3 (k = 2, 3, ..., 6) oraz p = 4 (k = 4, 5, 6) prze- działy absolutnej stabilności metod RKR są krótsze od odpowied- nich przedziałów absolutnej stabilności metody RK realizowanej z krokiem h/k,

2° dla p = 2 (k = 2, 3, ..., 6) oraz p « 4 (k = 2, 3) przedziały absolutnej stabilności metod RKR są dłuższe od odpowiednich prze-

działów absolutnej stabilności metod RK realizowanych z krokiem

bilności metody RKR jest jednostajny ze względu na k (k = 2, 3,

Dla metod RK rzędu p i p > 4 długość przedziału absolutnej sta- bilności zależy od parametrów metody# Aby badać stabilność metod RKR dla p > 4, trzeba rozpatrywać konkretną metodę RK.

5. RZ£D METOD JEDNOKROKOWYCH PO EKSTRAPOLACJI RICHARDSONA

W paragrafie 3 rozważano metody jednokrokowe Bobkowa rzędu p i m- -stopniowe z m > p. Metody te skonstruowano tak, aby miały możliwie najdłuższy przedział absolutnej stabilności. Powyżej zbadano zacho- wanie się przedziałów absolutnej stabilności metod RK rzędu p = 1, 2, 3 i 4 przy p = m, po zastosowaniu ekstrapolacji Richardsona.

Obecnie zostanie przedstawiony taki sposób doboru parametrów meto- dy jednokrokowej <1 z m>p, aby metoda <£R (metoda $ po ekstrapo- lacji Richardsona) miała maksymalny rząd dokładności.

Na podstawie twierdzenia 1.2 otrzymuje się, że dla każdej meto**

dy jednokrokowej rzędu p wielomian w(z) spełnia związek h/k,

3° dla p = 1, 2, 3 i 4 przyrost długości przedziału absolutnej sta-

(5.1) w(z)-ez = 0(zp+1).

Podobnie wielomian odpowiadający metodzie <£ użytej dwukrot- nie z krokiem h/2, spełnia związek

(5 .2 ) w ^ - e 2 = 0 (ZP+1) ,

ponieważ rząd metody tak realizowanej jest nadal p.

Po zastosowaniu ekstrapolacji Richardsona, wielomian w^Cz), od- powiadający metodzie $R i otrzymany na podstawie wzoru (4,3) z k =

=> 2, spełnia*związek ,

(5.3) w *(z)-ez = 0(z^+^ ),

Wynika to stąd, że metoda $ R jest rzędu p+1, gdy metoda jest

rzędu p. *

Poprzez odpowiedni dobór parametrów metody $ , można otrzymać w wielomianie w (z) określonym w (1.6) wartości d^ = 1 dla i = p+1, p+2, ..., m. Przy takich wartościach współczynników d^, dla wielo- mianu w(z) zamiast (5*1), (5.2), zachodzi odpowiednio

Tak skonstruowana metoda 5 jest oczywiście rzędu p. Wielomian sta- bilności odpowiadający tej metodzie ma rząd zgodności m. Metoda <£

jest rzędu m w przypadku użycia jej do rozwiązywania równań linio- wych.

Ekstrapolacja Richardsona jest dokonywana dla metody rzędu p, dlatego przy d^ = 1 otrzynuje się związek (5.3), zamiast związku postaci

w(z)-ez = o(zm+"*) ,

(5.4) w*(z)-ez = 0(zm+2).

Wynika to stąd, że ekstrapolacja Richardsona podnosi rząd metody o 1, zatem wielomian w^Cz) ma postać

-i z^ SP ,x z1

vTCz) - 1 + Z + . . . + p , + (p + 1 ) i + Z di ^{I T »} i=p+2

gdzie dp+2 £ 1. Aby dla wielomianu w(z) z dj, = 1 (i * p+1, p+2, ..., m) otrzymać (5.4), należałoby wykonać ekstrapolację Richardsona, przyjmując rząd równy m, to z kolei w przypadku dowolnego zagadnie- nia (1.1) nie podniesie rzędu dokładności rozwiązania.

Nasuwa się idea współczynników d^ tak, aby dla wielomianu w (z) określonego dla danej metody zachodził związek (5.4). Metoda jednokrokowa 5 skonstruowana tak, że odpowiadający jej wielomian spełnia (5.4), daje rozwiązanie z błędem rzędu 0(hm+^) dla za- gadnień liniowych.

Na podstawie zależności (2.1), jakie zachodzą pomiędzy stopniem metody RK a maksymalnym rzędem metody, naturalne jest dobrać d^

(i = p+1) dla p * 4, 5, 6 oraz di (i = p+1, p+2) dla p * 6, 7, 8.

Również w przypadku metod jednokrokowych RKB, gdzie na ogół jest m>p, można dobrać parametry metody tak, aby zachodził związek (5.4)*

Niech L (k = 2) oznacza przedział absolutnej stabilności meto- dy <£r.

Zachodzą następujące twierdzenia.

TWIERDZENIE 5.1. Dla metody 5 rzędu p = 1, 4, 5, 6 i (p+1)-stop- niowej trzeba przyjąć następujące wartości dp+.j» aby otrzy- mać związek (5.4). Metoda ta ma wskazane przedziały abso- lutnej stabilności L i L*:

Wartość p m

1 2

5 6

Wielomian w (z)

1+z+2 zf3 2!

i c _5 i a=0

Z

z\6 z

n 7 6T

i=0

(-3.0 0,0)

5 Z il+6 5T (-2.75,0)

(-3.25,0)

(-4.75,0;

(-4.25,0) (-6.37,0)

6 ^± y

6 7 Z ri+8 7T (-3.25,0) (-5.87,0)

□

TWIERDZENIE 5.2. Dla metody f rzędu p = 6, 7, 8 i (p+2")-stop- niowej, aby otrzymać związek (5.4) trzeba przyjąć d^+^, dp+²» które spełniają następujące zależności. Metoda ta ma

Wartość p m

wskazane przedziały absolutnej stabilności L i L : Związki pomiędzy L

8

7 9

8 10

^p+1 * %+2 -dg+8dy = 7 3dQ+12dy = 14

-dg+9dQ « 8 10dg+45d8 » 52

-diQ+10d9 = 9 11d10+55dg = 63

K

(-4.75,0) (-7.50,0)

(-5.12,0) (-8.25,0)

(-5.37,0) ^ (-9.00,0)

□

Dowód powyższych twierdzeń polega na sprawdzeniu, że przy tak ustalonych wartościach d^+>j oraz d^+^ i otrzymuje się związek

(5.4). Przedziały L i L* wyznaczono za pomocą maszyny cyfrowej.

PRZYKŁAD 5.1. Dane Jest zagadnienie początkowe

y' = -*y» yCo) = 1,

którego rozwiązaniem dokładnym Jest funkcja

Zagadnienie to rozwiązano metodą 5 rzędu z przykładu 2.1 i wyzna- czono rozwiązania w punkcie x * 1.0 przy wskazanych wartościach 0 ([¹²]).

6-

A £ =10-9

W

-6 1.1010-9 628 -2.9810-10 509

-1 1.3810-10 118 -1.97io-11 101

1 -1.9210-10 118 -4.2810-11 101

6 • -1.2010-9 610 -2.4010-10 525

Powyżej podano błędy względne (zadana dokładność £ * 10"^ oraz liczbę obliczeń wartości funkcji f(x,y). Dla ćJ = zachodzi w*(z)-ez = 0 (z8) . Przy metoda z przykładu 2.1 charaktery- zuje się długim przedziałem absolutnej stabilności L = (-5.62,0)

6. METODY JEDNOKROKOWE PIERWSZEGO RZ§DU ZORIENTOWANE NA ZAGADNIENIE POCZĄTKOWE

Rozpoznanie sytuacji w trakcie obliczeń na maszynie cyfrowej, kiedy wielkość kroku całkowania jest zdeterminowana długością przedziału absolutnej stabilności metody $ , a kiedy, zadaną dokładnością (rys.

3.1), jest ważne ze względów praktycznych. Otrzymanie informacji o rodzaju zależności kroku całkowania pozwala np. „przełączyć" algo- rytm na metodę jednokrokową o lepszych własnościach absolutnej sta- bilności (o dłuższym przedziale absolutnej stabilności), niż ma me- toda ?, a przez to zmniejszyć koszt otrzymania rozwiązania zagadnie- nia początkowego (1.1).

Tutaj jest przedstawione jedno z możliwych rozwiązań automatycz- nego wykrywania charakteru zależności kroku całkowania. Wykorzystu-

je się w tym celu dwie metody jednokrokowe pierwszego rzędu 2-stop- niowe :

Metoda należy do klasy netod (3.2), metoda ^ należy do klasy me- tod (3.3). Dla metodywielomian stabilności ma postać

W( Z) s 1 + Z + —z2 •

Przedział absolutnej stabilności wynosi (-3.00,0). Metodzie^ Opo- wiada wielomian stabilności

w (z) = 1+z + ⁵2 ,

przedział absolutnej stabilności wynosi (-6.00,0). Dla metody ^ zachodzi (5.4), tj.

A z ) -ez « 0 (z^) ,

natomiast dla metody lig zachodzi (5.3), tj.

w*(z)-ez * O(z^).

Metoda ta ma natomiast dwa razy dłuższy przedział absolutnej sta- bilności, niż metoda il^.

Jako podstawy w obliczeniach używa się metody 1-j. Metodą wyz- nacza się dwa rozwiązania: rozwiązanie obliczone z krokiem h oraz rozwiązanie obliczone z krokiem h/2. Wykorzystując wartości funkcji

f (x,y) obliczane przy realizacji metody^, można „skonstruować"

(bez obliczania wartości funkcji f (x,y)) również dwa rozwiązania me- todą 3ec*no 2 krokiem h i drugie z krokiem h/2.

Na podstawie rozwiązań otrzymanych m e t o d ą i metodą ^2 wyba- czane są dwie wielkości kroków całkowania. Za rozwiązanie przyjmuje się wartość (uzyskaną po zastosowaniu ekstrapolacji Richardsona) o- trzymaną tą metodą jednokrokową, która dopuszcza użycie większego kroku całkowania. Koszt równoczesnego stosowania metodyIł-j i §ł2 ^7“

nosi 5 obliczeń wartości funkcji f(x,y) na przedziale [x,x+h] .

PRZYKŁAD 6.1 ([12]). Zagadnienie początkowe

y' * 10 sgnsin (20x)y2, y2(°) = °#

y' - -10 sgnsin(20x)y1, y1(0) - 1,

ma rozwiązanie dokładne

y>j Cx) * | sin 10x|, y2(x) * I cos 1°x I •

Poniżej podano błędy względne ((y-y(x))/y(x)) otrzymanych rozwiązań w punkcie x = 1.0 oraz liczbę obliczeń wartości funkcji f(x,y) (ozna- czenie: [fj) wykonanych w celu otrzymania wartości w x = 1.0. Obli- czenia wykonano z zadaną dokładnością £* Dol:)<^r kroku całko- wania w metodzie 3^, i&2 oraz metodzie (realizowanej według powyższego opisu) przebiegał jak w pracy [²] .

Metoda 5 ^ Metoda 2

x _W W

1.0 -7.6410-4 3978 -6.6610-4 3346

Metoda

x

W

1.0 1.7910-3 1598 -2.7010-4

□

PRACE CYTOWANE

[ij V.V. Bobkov, Ob odnom sposobe postroenija odnoSagovych pravll pribliSennogo reSenija differencijal'nych uravnenij, Izv. AN BSSR, ser. fiz.-matem. No 4(1967.) (praca cytowana za [5]).

[2] J.S. Chomicz, A. Olejniczak, M. Szyszkowicz, Wyznaczanie kroku całkowania dla metod jednokrokowych i kwadratur, Raport Nr N-35 Instytut Informatyki Uniwersytetu Wrocławskiego, Wrocław 1978.

[3] P.J. Van der Houwen, Construction of integration formulas for initial value problems, North Holland, Amsterdam 1977.

[4] R. Jeltsch, D. Nevanlinna, Stability of explicit time discreti- zations for solving initial value problems, Report No 30, Feb.

1979, Depart, of Math. Faculty of Technology, University of Ou- lu, Finland.

[⁵] V.I. Krylov, V.V. Bobkov, P.I. Monastyrnyj, Vygislitel*nye me- tody, t. 2, Nauka, Moskva 1977.

[6] -, Vyćislitel,nye metody vysSej matematyki, t. 2, Vy5ej5aja Skola, Minsk 1975.

[7] J.D. Lambert, Computational methods in ordinary differential equations, Wiley, London 1973.

[8] J.D. Lawson, An order five Runge-Kutta process with extended region of stability, SIAM J. Numer. Anal. 4(1966), 593-597.

[9] -, An order six Runge-Kutta process with extended region of stability, SIAM J. Numer. Anal. 4(1967), 620-625.

[¹⁰] H.J. Stetter, Analysis of discretization methods for ordinary differential equations, Springer-Verlag, Berlin, Heidelberg, New York 1973 (tłum. ros. „Mir11 Moskva 1978).

[11] M. Szyszkowicz, 0 stabilności pewnych metod jednokrokowych, Raport Nr N-94, Instytut Informatyki Uniwersytetu Wrocławskie- go, Wrocław 1981.

[12] -, Metody jednokrokowe o podwyższonym rzędzie dokładności. Ra- port Nr N-118, Instytut Informatyki Uniwersytetu Wrocławskiego, Wrocław 1983.

[13] Wyznaczanie kroku całkowania dla metod jednokrokowych Bob- kowa. Mat. Stos. XXV(1984), 129-172.

[14] W. Świątek, T. Świątek, Automatyzacja obliczeń analitycznych na wielomianach wielu zmiennych, Instytut Informatyki Uniwer-

sytetu Wrocławskiego, praca magisterska, Wrocław 1981.

[15] Modem Numerical Methods for Ordinary Differential Equations (edited by G. Hall and J.M. Watt), Oxford 1976.