Wyznaczanie kroku całkowania w metodachjednokrokowych Bobkowa*

(1)

RO CZNIK I POLSKIEGO TO W ARZY STW A M ATEM A TY CZN EG O Seria III: M ATEM A TY K A STOSOW ANA XXV (1984)

M I E C Z Y S Ł A W S Z Y S Z K O W IC Z

Wrocław

Wyznaczanie kroku całkowania w metodach jednokrokowych Bobkowa*

(Praca wpłynęła do Redakcji 18.10.1982)

i. wsTęp

V niniejszej pracy rozważa się następujące zagadnienie począt- kowe

(1) y / s f (x,y(x)), x € [a,bJ, (2) y(a ) = ya ,

gdzie

y (x ) ,ya , f £ RS , ^{f tSi— *} RS , £2 = [a,b]XR , s>1 • Zakłada się, że

(3 ) fec(^),

* MR 1.1/04

[129]

(2)

130 M.SZYSZKOWICZ (b) V / \ Ilf (x,y)-f (x,y*)|| ^ L II y-y* || .

L>0 (x,y), (x,y*)eti

Przyjęcie założeń (3), ( ) zapewnia istnienie i jednoznaczność rozwiązania zadania (1 ), ( 2 ). Rozwiązanie to jest oznaczone przez y(X).

¥ pracy niniejszej rozważa się jawne metody jednokrokowe rozwiązywania zagadnienia początkowego ( 1 ), ( 2 )• Metoda jedno- krokowa jest najczęściej określona przez swoją funkcję przy- rostu $ f (x,y,h/), gdzie mając xn , yn , hn , kolejne yn+1 » przy- bliżenie wartości y(xn+hn), otrzymuje się na podstawie wzoru

\

( 5 } yn+1 = yn+hn (xn 'yn ’V •

Wielkość hn jest nazywana krokiem całkowania metody jednokro- kowe j. ¥ dalszej części tej pracy funkcja $^(x,y,h) jest ozna- czana przez i metoda jednokrokowa jest utożsamiana z funk- cją ^ • Wartości {yn } otrzymane z metody $ są nazywane wartoś- ciami rozwiązania przybliżonego. Dla podkreślenia zależności wartości yn od kroku h, przyjmuje się również oznaczenie y(xn ,hn_.j ). Dalej w pracy, wskaźnik n przy x oraz h będzie o- puszczany.

¥ trakcie rozwiązywania zagadnienia początkowego metodą $ dobierany jest krok całkowania. Dana metoda doboru kroku cał- kowania stosowana przy rozwiązywaniu zagadnienia ( 1 ), (2 ) wy- znacza siatkę

Djj = ■|x^€Rt xi+-| = xo = a » ^ = ^»***» ^

(3)

WYZNACZANIE KROKU CAŁKOWANIA W METODACH

JEDNOKROKOWYCH BOBKOWA 131

Wartości (yn }^_o rozwi^zania przybliżonego są określone tylko w węzłach siatki

W pracy zostanie przedstawiona metoda doboru kroku całko- wania , nazwana metodą podstawową doboru kroku całkowania, która to metoda może być stosowana dla dowolnych jawnych me- tod jednokrokowych. W metodzie tej wyznacza się dwa rozwiąza- nia przybliżone:

y(x+h,h) - wartość rozwiązania przybliżonego w punk- cie x+h otrzymana ze stosowania metody $ z krokiem całkowania h,

y(x+h,h/k) - wartość rozwiązania przybliżonego w punk- cie x+h otrzymana z k-krotnego stosowa- nia metody $ z krokiem całkowania h/k (k = 2, 3,..., 6).

Na podstawie tych wartości wyznacza się krok całkowania do na- stępnego etapu obliczeń przy rozwiązywaniu zagadnienia (1 ),

( 2).

W przypadku gdy metoda $ wymaga 1 obliczeń wartości funkcji f(x,y), to koszt realizacji metody podstawowej doboru kroku całkowania (w jednym etapie) wynosi (k+1)l obliczeń wartości funkcji f(x,y). Koszt ten jest zatem duży. Nie rezy- gnując z samej metody doboru kroku całkowania, można wykorzys- tać własności metody jednokrokowej $ ( w tej pracy są to me- tody Bobkowa), aby koszt ten obniżyć.

W § 2 tej pracy przedstawiono metodę podstawową doboru kroku całkowania. Następnie podano metodę otrzymywania metod

(4)

132 M. SZYSZKO VIC Z

jednokrokowych Bobkowa. V § k przedstawiono cztery sposoby wy- korzystania własności metod jednokrokowych Bobkowa, które umoż- liwiają zmniejszenie kosztów realizacji metody podstawowej do- boru kroku całkowania. W § 5 przedstawiono wyniki z ekspery- mentów na maszynie cyfrowej.

2. METODA PODSTAWOWA DOBORU KROKU CAŁKOWANIA Wyprowadzenie metody podstawowej

Niech y^(x+h,h) oznacza wartość w punkcie x+h otrzymaną meto- dą $ z krokiem całkowania h oraz y^(x+h,h/k) oznacza wartość w punkcie x+h otrzymaną poprzez k-krotne stosowanie metody z krokiem całkowania h/k, gdzie k może przyjmować jedną z war- tości {2,3,...}, natomiast i = 1, 2,.,,, s* Liczba s określa liczbę równań różniczkowych w zagadnieniu (1 ), ( 2 )• Do wyzna- czenia kroku całkowania h, który zapewnia otrzymanie wartości rozwiązania przybliżonego w punkcie x+h z żądaną dokładnością, zostaną użyte zależności między wartościami y^(x+h,h) oraz y± (x+h,h/k).

Dla metody jednokrokowej $ rzędu p zachodzi związek (6 ) y±(x+h)=yi (x+h,h)+C^°x)hP+1+c51^(x)hP+2+...,

gdzie y^(x+h) jest i-tą składową rozwiązania dokładnego za- gadnienia (1), ( 2 ) w punkcie x+h, (r)(x ) (r = 0, 1,...) są pewnymi funkcjami niezależnymi od h. Przyjmuje się, że

f £ C p+*Ca,b], a stąd cj^(x)cC^[a,b], Ponieważ interesujący jest jedynie rząd wielkości błędu, więc w dalszych rozważa-

(5)

WYZNACZANIE KROKU CAŁKOWANIA W METODACH 133 JEDNOKROKOWYCH BOBKOWA

niach we wzorze (6) uwzględnia się tylko składniki z najniż- szą potęgą h, tzn. z potęgą p+1• Znaku = dalej używa się do zaznaczenia, że zostały opuszczone wyrazy zawierające h w po- tęgach wyższych niż p+1•

Dla danej metody $ rzędu p i dowolnego punktu x zachodzi

Zależność (8) należy udowodnić.

Korzystając ze wzoru (7 ) z krokiem całkowania , otrzy- muje się

Stosując wzór ( 9) k —krotnie kolejno w punktach x, x+£ » •••t (7) yi (x+h)-y± (x+h,h ) = C± (x)hP+1

oraz

(8) yi (x+h)-y± (x+h,h/k ) = k ^ (x )(h/k )p+1 , gdzie C^(x ) =

(9)

h x+(k-l) dla danej metody $ dostaje się

i ...

(6)

134 M.SZYSZKOWICZ + £ $ f (x+(k- i ) ^ y T (x+( k - i ) | ) , | ) +ci (x+( k - D

%)[%f

,

gdzie y (x) = (y 1(x),y2(x;,...,yg (x);.T

Ponieważ opuszcza się składniki zawierające wyrazy z h w potędze wyższej od p+1, zachodzą związki

c±w i%r' = ci ( x+t ) ( i r ' = ••• = «s1(»(^ < > i)(ijp+1

Związki te dają

y± (x+h) i y± (x )+^ (x »yT(x)»i)+| M x+i»y T (x+E )»£)+

+•••+£ £f (x,yT (x+(k-l) ^)»|)+kCiCx )(t)

Zastępując wartości y± (x+^)» y±(x+^£ )» • • • » Yi(x+(k“ 1)

(i = 1, 2,..., s) odpowiednimi wartościami z (9) i zaniedbu- jąc wyrazy z potęgami h większymi od p+1 w rozwinięciach funkcji:

i / h T, h\ T, v / h ^p+1 h \

$f(X+k»y (x+k)+c <x >(ki »k b t ( 2h T/ 2h \ T( h\/h\p+1 h'l M x+"k»y (x+"k)+c (x+k A kI ’kh

§f (x+(k-l) £,yT (x+(k-l) |)+cT (x+(k-2)

(7)

ostatecznie otrzymuje się

yi(x+hj = yi ('x+h,|)+kCi(x)[|) ^,

co jest zgodne z (8),

Po odjęciu stronami zależności (8) od (7) dostaje się

„ , ₁

związek, który pozwoli wyznaczyć wielkości C^(x)hp , miano- wicie

(10) A ± = yi (x+h,h/k)-y± (x+h,h) = C± (x)hp+1 .

Przez krok H rozumie się taki krok całkowania dla danej meto- dy który zapewnia otrzymanie rozwiązania przybliżonego w punkcie x+H z żądanym błędem. Niech 6 oznacza wektor do- puszczalnych błędów. Żąda się, aby spełnione były nierówności

(11) |ci(x)HP + 1 | i | yi (x+H)-yi (x+H,H)|.

Kroku całkowania H szuka się w postaci H = ~, zatem należy wy- znaczyć wielkość Ca). Krok h jest znany i jest to krok całkowa->

nia aktualnie stosowany w obliczeniach. W (11) zakłada się, że główny składnik błędu dobrze przybliża cały błąd. Ponieważ wyznaczone H nie zawsze to zapewnia, we wzorze (12) "poprawia"

się H, uwzględniając niejako odrzucone składniki błędu.

Nierówność (11) można zapisać w postaci

(8)

136 M.SZYSZKOWICZ

£li.

| = | c i ( x ) ( g ) p+1|

i dalej na podstawie (10) dostaje się

|ci(*)hPł1 u U ± |-sL .

Stąd zaś

\A± \ kP

P > 1 ^K

i / . p+1 /.p UJ ( k - 1 ). \

Ostatecznie wzór na U) przyjmie postać

(

^{1 2}

)

CO = const

k^-1 max 1 £ i <£ s

Stała const (const>1) występująca we wzorze (12) zapewnia poprawne wyznaczanie kroku H, Przyjmuje się H mniejsze niż to wynika z (11).

Błąd rozwiązania otrzymanego z krokiem całkowania H jest rzędu ^o (^h^+M . V trakcie wyprowadzania wzoru (12) jako wartość początkową przyjmuje się dokładną wartość rozwiązania zagad- nienia początkowego, tj. wartość, której rząd dokładności jest co najmniej o (h^+^) • ¥ celu otrzymania takiej wartości do na- stępnych etapów całkowania, korzysta się z y^(x+H,H),

y^(x+H,H/k) i po zastosowaniu ekstrapolacji Richardsona okreś- la się wartość

(9)

y ...y^(x+H,H/k)-y^(x+H,H)

WYZNACZANIE KROKU CAŁKOWANIA W METODACH 137 JEDNOKROKOWYCH BOBKOWA

y*(x+H,H ) = yi(x+H,H/k ) +

k - 1

Otrzymana wartość y^(x+H,H) jest obarczona błędem 0(H^+^) i może być zatem przyjęta jako wartość początkowa do następnego etapu obliczeń.

Takie postępowanie pozwala otrzymać na końcu przedziału całkowania wartość, która w zasadzie nie różni się od rozwią- zania dokładnego o zadany z góry wektor błędów 6 . Metoda do- boru kroku całkowania, tutaj przedstawiona, została nazwana metodą podstawową. Inne metody podane w tej pracy są modyfika- cjami powyższej metody i głównie chodzi w nich o zmniejszenie kosztów otrzymywania wartości y^(x+h,h), y^(x+h,h/k). W dal- szym ciągu tej pracy w opisie innych metod nie będzie już wy- raźnie zaznaczone, że metoda stosuje się do układu równań róż- niczkowych, tj. opuszczony będzie wskaźnik i.

Schemat blokowy realizacji metody podstawowej doboru kroku całkowania

Algorytm, który realizuje przedstawioną powyżej metodę doboru kroku całkowania, może mieć postać jak na rysunku 1. W reali- zacji wzór (12) na liczbę W ma postać

OJ— 1.25

gdzie 6± = eps|y^(x+h,h)| oraz krok H dający żądaną dokład- ność eps jest równy H»k = ~ • Krok H jest tak dobierany, że całkowanie z krokiem H daje wartość rozwiązania przybliżonego

l * i l

k(k^-i)eps 1-£i-£s ly^fx+hjh)max

(10)

138 M.SZYSZKOWICZ

(11)

z błędem większym niż dopuszczalny, a całkowanie z krokiem H/k daje wartość rozwiązania przybliżonego z błędem mniejszym niż dopuszczalny* W tym sensie krok całkowania H może być uwa- żany za optymalny. Wykonanie ekstrapolacji na otrzymanych wartościach rozwiązania daje dokładniejsze rozwiązanie. Zatem dla procedur opracowanych według powyższego schematu zachodzi proporcjonalność błędu globalnego do wielkości eps, co można zaobserwować na wynikach obliczeń prezentowanych w § 5 tej pracy. Metoda podstawowa doboru kroku całkowania z k = 2 była publikowana w [1].

Najczęściej w realizacji powyższej metody doboru kroku całkowania przyjmuje się k = 2, co w przypadku metody $, któ- ra potrzebuje 1 obliczeń wartości funkcji f(x,y) w jednym kro- ku, wymaga obliczenia 31 wartości funkcji f(x,y) w celu wyzna- czenia liczby CO. Dla dowolnego k potrzeba (k+l)l obliczeń wartości funkcji f(x,y). Dalej zostaną przedstawione sposoby zmniejszania kosztów otrzymania liczby CO, w przypadku stoso- wania metod jednokrokowych Bobkowa.

3. METODY JEDNOKROKOWE BOBKOWA

Bobkow w pracy [ 2 ] przedstawił następujący sposób otrzymania metod jednokrokowych.

Całkuje się równanie (1 ) na przedziale [x,x+h] i dostaje się równość

y(x+h) = y(x) + x+h

j

f(t,y(t))dt.

WYZNACZANIE KROKU CAŁKOWANIA W METODACH 139 JEDNOKROKOWYCH BOBKOWA_____________________________________

(12)

1^0 M.SZYSZKOVICZ Po dokonaniu zamiany zmiennych t = x+GCh przyjmuje ona postać (13) y(x+h) = y(x)+h

Do obliczenia wartości całki w (13) używa się kwadratury, która po wprowadzeniu oznaczenia z(a.) = f (x+a. h,y(x+01 h )J jest postaci

I z (a) da = ^{^ 2} A. z (ot.)+r ,

0 i=0

gdzie r jest błędem stosowanej kwadratury.

Można teraz (13) zapisać w formie równości

(14) y (x+h ) = y(x;+h ^ A.z ( Ot, )+r • i= 0 1

Wyboru parametrów ^ = ) dokonuje się

tak, aby zastosowana kwadratura miała możliwie wysoki rząd do- kładności. Rezygnuje się jednak z najwyższego rzędu dokładnoś- ci w celu uzyskania swobody w doborze parametrów. W tej pracy żąda się, aby zastosowana kwadratura była dokładna dla wszyst- kich możliwych wielomianów stopnia co najwyżej p-1, gdzie

0<p^2q+2 i p jest założonym rzędem metody jednokrokowej. Pro- wadzi to do następującego układu p równań z 2q+2 niewiadomymi {_ } t (i = 0» 11 • • • t q);

J f (x+ ^CLh,y (x+ ^ah)) d a •

(15) £2- ś.^=1, ^ ^i^i3 i+1 (J=1 t 2,..., p-1).

i=0 i=0 J+

(13)

WYZNACZANIE KROKU CAŁKOWANIA W METODACH JEDNOKROKOWYCH BOBKOWA

Wiadomo z teorii numerycznego całkowania (np. § 5.4.1 i

§ 5.4.2 w [2]), że kwadratura zastosowana w (14), mająca naj- wyższą algebraiczną dokładność 2q+1, może być utworzona i to

jednoznacznie dla dowolnego q£.0. Przy p = 2q+2 układ (15) nia jednoznaczne rozwiązanie.

Jest to rozwiązanie, jakie otrzymuje się przy wyznaczaniu kwadratury Gaussa. 0 ^ - zera (q+1)-szego wielomianu ortogonal- nego w przedziale (0,1 ) z wagą 1. A ,> 0 oraz ^ A. = 1} stąd

1 i=0

0< Ai ^ 1. Natomiast dla 1« p < 2q+2 układ (15 ) ma przynajmniej jedno rozwiązanie, czyli można skonstruować metodę jednokro- kową Bobkowa dowolnego rzędu.

Dla oszacowania błędu r w (14) łatwo zauważyć, że układ (15) można otrzymać z rozwinięcia względem potęg h wyrażenia

y(x+h)-y(x) oraz

i h & A.y ( x+ Ot h )

i=0

tak, aby rozwinięcia te były zgodne włącznie do wyrazów zawie*

rających h^. Taki sposób otrzymania układu (15) pozwala napi- sać następującą postać błędu r w (14): *

(,6 ; Ł ).

Tak więc, gdyby wszystkie wartości y(x+0(^h) (i=0, 1,..., q) znane były dokładnie, wtedy z (l4) można byłoby otrzymać roz- wiązanie y(x+h) z błędem rzędu o(hp+^j postaci (16). Wartości

(14)

M.SZYSZKOWICZ

y(x+ OŁ^h) można wyznaczyć w sposób przybliżony i, jak widać z (14), wystarozy wyznaczyć je z błędem rzędu o(hp ). Zamienia- jąc h na OŁ^h, można podobnie do (14) napisać wyrażenie dla y(x+ a ±h)

y(x+ o6±h) s y(x;+h *i Bjz (fy a i)+ri » j—0

gdzie z ( ) s f (x+ jb Otih,y(x+jb oc^h)). ¥ celu znajdowania wartości y(x+ Jbj O^h) można znów skorzystać z analogicznej równości przybliżonej jak w (14). Mając na uwadze, że przed sumą w (14) stoi czynnik h, każde następne przybliżenie może być wyznaczone z mniejszą dokładnością. Tak więc, wartości y(x+ OLjh.) wyznacza się z lokalnym błędem rzędu 0(hP ), wartości y(x+CX^ Już tylko z lokalnym błędem rzędu o(hP” *) itd.

Jest to równoznaczne z odrzuceniem na każdym kolejnym etapie ostatniego równania z układu (15 )• ¥ związku z obniżeniem do- kładności zmniejsza się również liczba parametrów stosowanej kwadratury. Ostatnimi równościami, na których proces się za- trzymuje, są

y(x+ OŁ± Pj... ~ y ( x ) + a ± Jbj... y 1hf(x,y) .

Błąd tych przybliżonych równości jest rzędu 0(h ) i są to wzo-2 ry prostej kwadratury lewych prostokątów.

¥ trakcie realizacji tych metod jednokrokowych proces przebiega odwrotnie. Najpierw wzorami określającymi metodę pierwszego rzędu wyznacza się wartości w odpowiednich węzłach.

(15)

Korzystając z uzyskanych wartości,konstruuje się rozwiązanie o rząd wyższej dokładności i tak aż do wzorów dających wartość przybliżoną dla y(x+h) z błędem rzędu o(h^+^). Bobkow określa swoje metody jako metody jednokrokowe rekurencyjnego podwyż- szania rzędu.

Algorytm uzyskiwania metod Bobkowa

W metodach jednokrokowych Bobkowa) w celu obliczenia wartości y(x+h), stosuje się kwadraturę (metodę numerycznego całkowa- nia) postaci (14)

q ,

y(x+h) = y(x)+h A.y (x+ a.h)+r, i=0

jdzie r jest błędem kwadratury. Dla metody rzędu p, parame- try {A± } , muszą spełniać następujące równania

A . = 1, OŁ. = k+T = 2,.,., p-1;

i=0 1 i=0 1 1

p ^ 2ą+2)•

Często parametry {A^| nazywane są wagami, a parametry węzłami kwadratury.

Ponieważ wartości y(x+0(ih) występujące po prawej stronie w (14) nie są znane, należy je wyznaczyć. Można to zrobić za pomocą kwadratury, która jest zastosowana w (14). Narzuca się pewien proces rekurencyjnego otrzymywania kolejnyoh wartości występujących w (14). Proces ten daje się przedstawić w posta-

(16)

ci algorytmu, który realizuje się etapami.

Oznaczenia:

1 M. SZYSZKOVICZ

- węzeł o numerze i w etapie r, A^ (r) - waga o numerze i w etapie r,

*ri - liczba węzłów (wag) użyta do wyznaczenia war- ( X*— 1 )

tości w węźle OL: w etapie r.

Etap 1, Na podstawie (14) i (15) otrzymuje się wartość rn+ 1

y = y +h

J n+1 J n i=0 ‘— i. i J n+ai^ (1 )

oraz warunki

r ' a (1) i ^ 1

c Ai = 1, C A± (0Ł± ) =

i=0 i=0

(k = 1, 2,,.., p-1; p^2q+2).

Spełnienie powyższych równań dla |a^ ^ | oraz gwarantuje uzyskanie wartości yn+^ z błędem rzędu o(h^+^).

Etap 2, Wyznacza się wartości y ze wzorów n+(Xi

(17)

q2i

= yn+«i,)h c *!2,y

J = 0 n W 1 1 ^ 2 ’

Ł » '-* '1 ^ 4

1 i J

Aby otrzymać y ^ ^ z błędem rzędu O(h^), muszą być spełnio- n+a^

ne następujące równania:

q2i q2i

E * f = , . E A<2,( a (2> f = ^

i=6 J j=o 11+1

Z równań tych otrzymuje się wartości odpowiednich węzłów i wag.

Etap r. Wyznacza się wartości y (r-1 ) ze wzo- n+“i a j •••a o

rów

y (1) (2) (r-1) = y (1) (2) (r-1 ) x n+01^ CXj * • • n+Ol^ Otj • • •

v w a

* h ^ s y (1) (2) (r-1) (r) • n+a± ad ...a0 0Łfl

Aby otrzymać y (1 ) (2 ) (r-1) z bl9dem rzędu 0(hp+2*"r ), n+a, _i a . • • • _j _o

muszą być spełnione następujące równania:

(18)

A (r) - 1 *¥* A Cr)('OLCr) F - -1- t lS s " ’ s C s J - k+1

S =U S=U

(k s 1, 2,..., p-r; p-r+1 Ś 2qrQ+ 2 ).

Podobnie jak w poprzednich etapach z równań tych otrzymuje się wartości odpowiednich węzłów i wag.

Etap p (ostatni). Wyznacza się wartości y (1) (2) (p 1) n+a. _i a. ...a _j _u

ze wzorów

(1) (2) (p-1) n+OC^ ^j * * *^u

qpu - n W 1M 2)— ^r ^{1 ) *} v=0

Wzory te są uzasadnione tym, że stosuje się przybliżenie war—

tości y £,,) (p-1) z rzędu 0(h ). Do uzyskania2 nł0li “ j ---“u

takiej dokładności rozwiązania wystarcza metoda Eulera* Speł- nione muszą być równania

’£ ^ = 1 , , = 0 .

^ V ’ qpu

To kończy proces wyznaczania węzłów {&} i wag {a |i otrzymu- je się wzory określające metodę rzędu p.

(19)

JEDNOKROKOWYCH BOBKOWA 1 k7

W trakcie realizacji metod jednokrokowych Bobkowa, powyż- szy algorytm jest realizowany od końca. Zaczyna się więc od etapu p, w którym oblicza się wartości z błędem rzędu 0(h ). 2 Kończy się na etapie 1, gdzie dzięki kolejnemu podwyższaniu (na etapach p, p-1,..., 2) rzędu dokładności wyznacza się war- tość yQ+^ z błędem rzędu 0(h^+^).

Układy równań na wagi {a } oraz węzły {&], otrzymywane na każdym etapie powyższego algorytmu, przy przyjętych ogranicze- niach na rząd rozwiązania oraz liczbę węzłów i wag, mają zaw- sze przynajmniej jedno rozwiązanie.

W praktyce celowe jest zmniejszenie liczby węzłów i wag w kolejnych etapach algorytmu. Uzyskuje się przez to obniże- nie stopnia metody Bobkowa. Z drugiej strony jednak większa od wymaganej liczba węzłów i wag pozwala swobodniej dobierać ich wartości.

Układy równań na {a } i {oc} otrzymywane w kolejnych eta- pach algorytmu, można wyprowadzić w inny sposób, przedstawio- ny poniżej.

Rozwija się wyrażenie L = y(x+h)-y(x) oraz wyrażenie R s h ^ A.^ \y‘ (x + a / ^ h ) względem potęg h i żąda się zgod-

i=0

ności tych rozwinięć do wyrazów zawierających h^.

Dla L otrzymuje się rozwinięcie L = hyf ( x ) + y ;/(x)+ ••• +

(20)

I

1*ł8 M.SZYSZKOWICZ

'■ l'?fud'ł^J' o^wauu, Litx y f / ^ ^ wstawia się wyrażenia

Wyrażenie R rozwija się w następujący sposób. Za y'( x+a±( h)

y'(x+cJ1^h ) s y ł( x ) +

+ oc^i;h Y 2 Aj2)y w(x+a(1)a (2)h)+0(hp ) ^.

j=0 3 3

Po dokonaniu tego wstawienia wyrażenie R przyjmie postać

R = h C A ^ i y (x ) + i=0

<łoi

+ A?2)y»(x+aC1 )a !2)h ) + 0(hP ) J .

j_0 3 3 J

Następnie za y M(x+CX^1 ^aj2>h ) wstawia się

y"( x+tt^1 *(xj2>h ) = y"(x ) +

a<3V " ( « a < 1>a<2,a < 3)h) + o(hP"1 ).

3 1=0 3 L

Teraz wyrażenie R jest postaci

R = h ^ A^l)[y'( x)+a±(1)h Af2)[y"(x) +

i=0 X 1 j=0 J

A^3)y "/(x+«i(1)otj(2)cy10 ) h)+o(hp~1)] +o(hp)].

(21)

I WYZNACZANIE KROKU CAŁKOWANIA W METODACH 1^9

JEDNOKROKOWYCH BOBKOWA

Wskazane powyżej podstawienia wykonuje się p-1 razy. Po p-1 krokach wyrażenie R jest następującej postaci:

q (1 )r , R = h T2 A, [y ( x) +

i=0

j=0 J 1=0

\v=0

+ 0(h)] + o(h2)] . . . + 0(hP ) ] . p-1

Porządkując w R wyrazy względem potęg h i przyrównując L do R ? otrzymuje się warunki na {A} i {&}

i=0 * * j=0 J

£ - *1 .

i=0 1 .1=0 J 1 J 1=0 x

i=0£ A (1Vi i1}

q2i

j=0z : A<2 V (1>a (.2) J ^ J ^q£₁₌₀ (3)

(22)

a (1V 2)Ct(3> QL(r“ 1)- V'* A (r)____ 2- ...a± aj ...a0 Jz© » “ r!

& C £ A<3 >...

i=0 1 1 j=0 J J 1=0

(1) C2) C3) ^(P-1). A (P) _ _i . ...cy± a , ...a u > Av - Pi •

Po wyłączeniu i pogrupowaniu węzłów -[oc ]■ przed sumami, otrzy- muje się następujący układ (właściwie jest to układ złożony z równań typu (15 ))i

£ aa(1) = 1, i=0

q\

j q

£ A f •1 = 5 f , i=0 1 .1=0 J

£ c A f > < a < 2 ) r 2 * « > . . .

i=0 j=0 J J 1=0

q ( r - 11n r _ i ) ( r )

^ n 0<o s=0 As - rl ’ o =0

(23)

£ £ A f > < a f V - 2

1=0 j=0 J J 1=0

Q[CP-1) . V A(p) - -*• .

• • u ^ V “ p!

Łatwo zauważyć, żo ten układ równań spełniają parametry {et}, { A j otrzymane w kolejnych etapach algorytmu. Tak więc powyższy

układ ma rozwiązanie.

W szczególności| jeżeli do powyższego układu wstawi się rozwiązania na {cxj, {a } z 2, 3,..., p-tego etapu algorytmu,

to otrzymuje się układ (15)s

(24)

152 M. SZYSZKO WIC Z i=0t

1 1

(P-1) (P-2) Pl¹

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

yn+a - przybliżenie wartości y(xn+o(h),

k k

yn+ot “ przybliżenie wartości y(xn+0(h) z błędem 0(h ) i analogicznie

fn+a = f (xn*a h 'ya+a > oraz fn+ a =

PRZYKŁAD 1• Dla ilustracji poniżej jest wyprowadzona meto*

da jednokrolcowa Bobkowa drugiego rzędu* ¥ tyra wypadku układ (15 ) ma postać

£ a = i, £ a a = i.

i=0 1 1=0 1 1 d

Dla q = 0 otrzymuje się jednoznaczne rozwiązanie Aq = 1,

0(q = ^"2 -*■ prowadzi to do następujących wzorów określających me- todę drugiego rzędu:

(25)

WYZNACZANIE KROKU CAŁKOWANIA W METODACH 1/53 JEDNOKROKOWYCH BOBKOWA

Przy q = 1 układ (15) przyjmuje postać

A 0 + = 1 ’ ^ O ^ O + = H *

Wybierając OC0 s 0, ^ s 1, znajduje się, że AQ s A1 = ^ i otrzymuje się wzór

y 3 . = y3 + f (f3 + f3 J .

■’n+l Jn 2 ' n n+1 '

Jest to niejawny wzór trapezów.

Korzystając ze wzoru Eulera, można podać następującą jawną metodę drugiego rzędu opartą na powyższym wzorze :

y2 = y3 + kf3 »_łn+1 _n _n_f

y3 Jn+1 = y3 + | (f3 + f2 ).O Jn Z ' n n+17

Przy konstrukcji metod jednokrokowych Bobkowa jest duża dowolność stosowania wzorów dających rozwiązanie w węzłach {ot^j z żądaną dokładnością. Metody Bobkowa mają tę własność, >

że każdy wzór wchodzący w ich skład określa wartość rozwiąza- nia przybliżonego w pewnym punkcie przedziału [x,x+h]. Takiej własności nie mają np. metody Rungego-Kutty,’ gdzie wartości liczone w punktach pośrednich między x a x+h nie są na ogół wartościami rozwiązań zagadnienia początkowego w tych punk- tach.

(26)

M.SZYSZKOWICZ

k. DOBÓR kroku całkowania przy wykorzystaniu własnoSci METOD JEDNOKROKOWYCH BOBKOWA

Metoda pierwsza

Metoda jednokrokowa Bobkowa rzędu p jest określona przez ciąg wzorów, które pozwalają przez ich kolejną realizację podnieść rząd dokładności rozwiązań przybliżonych w węzłach {oC^}

(i = 0, 1,***, q). Ostatni wzór realizowany w metodzie Bobko- wa rzędu p zgodnie z 14 jest postaci

(17) ' yP+J = yp+1+h Jn+1 ■'n f \₁₌₀ A 4fp „ ,' i n+Ot, ₁ 1

gdzie A^, spełniają układ (15),

A

j

= 1 ■^‘iP^'i = "7+7 ^ =

i=0 i=0 J

Bardzo często w metodach jednokrokowych Bobkowa ustala się wartość węzłów OLq = O oraz = 1 . Pakt ten może być wyko- rzystany do zmniejszenia liczby obliczeń wartości funkcji

w realizacji metody podstawowej doboru kroku całkowa- nia* Zwykle dla wszystkich metod jednokrokowych, jeżeli w kon- strukcji metody $ występuje wartość fn , można wykorzystać tę wartość do otrzymania rozwiązań y(x+h,h) oraz y (x+h/k,h/k ), co zmniejsza ogólną liczbę obliczeń wartości funkcji f(x,y) do (k+l)l-1, Natomiast ustalenie CŁ^q = O oraz 0( = 1 dla metod jednokrokowych Bobkowa ma daleko idące konsekwencje*

Normalnie, tzn* tak jak się oblicza w przypadku dowolnej metody jednokrokowej przy wyznaczaniu wartości y(x+h,h)

(27)

JEDNOKROKOWYCH BOBKOWA______________ 155 zamiast wzoru (17) realizowany jest wzór

<«> C ! = C 2 ♦ * ( v T a + c v L , ) •

gdzie wartość y^+2 została otrzymana na podstawie ekstrapo- lacji Richardsona. Taki sam wzór realizowany jest przy wyzna- czaniu wartości y(x+h/k,h/k)91j•

Następne wartości y (x+2h/k,h/k ), y(x+3h/k,h/k )»•••»

y(x+(k— 1)h/k,h/k) oblicza się wzorem

(19) ^n+ (r+1)/k = yn+r/k + k (A0fn+r/k+JŁ ‘k±fn+ (r+o^) /*) (r : 1, 2, • • • | k— 1 )•

We wzorach (18) i (19) występuje wartość f^+2» fn+ 1 ’ Podczas gdy dla otrzymania wartości y^*j wystarczy wartość fP , tak jak we wzorze (17)* Zatem rząd metody Bobkowa nie zmieni się, jeżeli przy wyznaczaniu wartości y(x+h,h) oraz y(x+h/k,h/k) zamiast (18) zastosuje się wzór

(20) yP+1 . s yp+2 + ^ A.fP /. (r = k, 1)

' y^+r/k n k i n+ro^/k v

(28)

156 M.SZYSZKOVICZ

oraz do wyznaczenia y(x+2h/k,h/k ), y (x+3h/k,h/k ), . . . , y (x+(k-1)h/k,h/k ) wzór postaci

(^{p i}) VP+1 _ _P+1 . h & A «p

yn+(r+1)/k “ ^n+r/k k i n+(r+oti)/k (r s 1, 2,. . . , k-1 ).

Użycie wzorów (20), (21) w realizacji metody podstawowej dobo- ru kroku całkowania pozwala zredukować cały koszt obliczeń do (k+1)(l-l) obliczeń wartości funkcji f(x,y).

PRZYKŁAD 2. Metoda Bobkowa drugiego rzędu postaci

2 3 3

y •'n+l 1 = y + hfJ ,■'n n 7

y3 - + — (f*3 + A n+1 n 2 ^ n n+l'»

realizowana w metodzie podstawowej doboru kroku całkowania z wzorami (20), (21) ma postać

2 3 2

y i = y + ńf ,n+1 J n n f

y^ = y^ + — (f^ + f^ ) □ yn+1 yn 2 ' n n+17, U

Metoda druga

Jak zostało pokazane w § k te’j pracy, błąd metody Bobkowa rzędu p może być przedstawiony w postaci (16):

(29)

r = - -1 t A ^ J + O O ^ ) .

W tyra przedstawieniu, jak i w trakcie wyprowadzania metod jednokrokowych Bobkowa, zakłada się, że znane jest rozwiąza- nie dokładne zagadnienia (1), (2) w punkcie x.

Można zażądać, aby parametry (l = °» 1»»*»» q) oprócz układu (15) spełniały dodatkowo

(22) (P+i)i ” pT ±& *1*1 = fl* »

gdzie £ jest odpowiednio dobraną, ustaloną liczbą rzeczywi- stą. W ten sposób zadana jest rodzina metod Bobkowa rzędu p, których błąd wyraża się w postaci

(23) r = hP+1J-y(p+1)(x)+0(hP+2).

Jeżeli ostatnim wzorem w ciągu wzorów określających metodę jednokrokową Bobkowa będzie wzór postaci (17)* to błąd nie będzie się wyrażał w postaci (23)* ponieważ wartość £ w głów- nym wyrazie błędu (23) zostanie zmieniona przez błędy wnoszo- ne wartościami y^+Qt (1 = 0, 1,..., q). Wystarczy zatem wyzna- czyć wartości y^ ^ (i = 0, 1,...,q) i zamiast wzoru (17) sto- n+o(^D+1 sować wzór

(24) vp+1 - vP+1 +h 5^ A fP+1yn+1 " yn i n+o^ *

(30)

aby błąd metody wyrażał się wzorem (23). V celu otrzymania wartości yf>+l (i = 0» n+u_^ q) można skorzystać z metod jed- nokrokowych Bobkowa lub innych metod jednokrokowych.

Dobierając odpowiednio parametry {0^ } (i = °»

•••i q)* można utworzyć metody jednokrokowe o różnych wartoś- ciach ^ w (23). Zakładając, że dana jest metoda jednokrokowa

, która w (23) ma wartość ^ można utworzyć metodę jed- nokrokową §2 , która z kolei w (23) ma wartość £ = ^ 0/kP (k = 2, 3,..., 6). Tak więc, wartość otrzymaną ze stosowania metody można traktować jako y(x+h,h), a wartość otrzymaną metodą $2 można traktować jako y(x+h,h/k).

PRZYKŁAD 3* Przyjmuje się rząd metody p = 2 i na podstawie (15) otrzymuje się układ

i=0 * i=0

oraz dodatkowo z (22)

^ = 5 - H '

Ustalając q = 1 oraz 0(Q = 0, dostaje się następujące rozwią- zanie :

(31)

JEDNOKROKOWYCH BOBKOWA 159

Z warunkami O < 0^ 1 dostaje się

- 75* i'<i ^•

Przyjmując ~ , otrzymuje się wzór postaci (24)A

(_{s J}^{2 5}) y3 < = y3 + § (f3 + _{J n+1} J n 2 ' n n+1 yi )*

Wartość 3 może być wyznaczona np. na podstawie wzorów

^ U

yn+l/2 “ yn + 2 n *

yn+1 = yn + ^ 0+1/2 #

Dla rozwiązania określonego w (25 ) zachodzi związek

y(x+h) i y3+1 - -Ą _{y (3)(x)h3}

i rozwiązanie y ^ jest interpretowane jako y(x+h,h). Chcąc3 otrzymać rozwiązanie yn+j» które można przyjąć za y(x+h,h/2), — O tzn. dla którego zachodzi związek

y (*♦*») = ? 3+1-a -Ą y <3>(x)(|) , 1 $O

należy przyjąć ^ a - j-g- = .

(32)

160 M. SZYSZKOWICZ Dla ^*= — ^ otrzymuje się wzór postaci (24)

<2Ć> yL i = ń ł 3 (fn + " L 3/4_{) *}

gdzie wartość o 3/^4 może **y^ wyznaczona np. za pomocą wzorów

y n + 3 / 8 = y n + f tfn -

V n * 3 A = y l * l “ n+3/8

Korzystając z wartości określonych w (25 )» (26), można reali- zować metodę podstawową doboru kroku całkowania i na podsta- wie y^ .. oraz y^ . określić wartość_J n+1 J n+1 ^

-3 _ 3 . -3 . yn+1-yn+1 n y (x+h,h ) = yn + 1 + ---- 3---- . U

¥ przykładzie 3 nie ma praktycznie oszczędności obliczeń wartości funkcji f(x,y) w porównaniu z wyznaczaniem w meto- dzie podstawowej wartości y(x+h,h) i y(x+h,h/2). Dopiero wy- znaczanie wartości » którą można przyjąć jako y(x+h,h/k) i k > 2, daje wyraźną oszczędność obliczeń funkcji f(x,y).

Żądanie, aby metoda dawała wartość z a meto- da $ g dawała wartość z , prowadzi do tego, że węz- ły {ot.^ (i = 0, 1,..., q ) dla ^ oraz $ 2 są na ogół różne*

Fakt ten nie pozwala wykorzystać wartości rozwiązań przybli- żonych wyznaczonych w do konstrukcji rozwiązań przybliżo- nych w metodzie $2#

(33)

Odrzucając warunek, że metodzie ^ odpowiada w (23) war- tość ^ a metodzie odpowiada w (23) wartość

^Q/kP , przyjmuje się, że metodzie ^ odpowiada w (23₎

^ ^ , a metodzie ^odpowiada w (23) = ^2* ^artośoi fi ^ (2 1 można ustalić, że | ^ | > Ifjjl*

Przy wyprowadzaniu metody ^ i 4^ można teraz zwrócić baczniejszą uwagę na węzły (1 = 0* <1)3- starać się tak je dobrać, aby było możliwie jak najwięcej wspólnyoh wartości rozwiązań przybliżonych w ^ i $,,•

Dla rozwiązania otrzymanego metodą ^ zachodzi związek

(27) y(x+h) = yjj£] + J'1y^P+1^ (x)hp+1

oraz dla rozwiązania otrzymanego metodą

(28) y(x+h) = yP*J + ^2y(P+1) ^x)hP+1

Postępując podobnie jak przy wyprowadzaniu metody podsta- wowej doboru kroku całkowania, można wykorzystać zależności

(27), (28) do wyznaczenia kroku całkowania, który zapewnia otrzymanie rozwiązania z żądaną dokładnością* ¥ tym przypadku wzór (12) na wielkość Có przyjmie postać

(34)

162 M. SZYSZKO VIC Z

gdzie const >1, A = yn+l” ^n+1* a wektor dopuszczalnych

b ł ę d ó w *

Na podstawie wartości yj^*j oraz yj^j konstruuje się roz- wiązanie przybliż me

D+1 -p+1

*/_ , , \ P+1 *n+1 ^n+1 y (x+h,h ) = y£+1 + — --- ,

£ l - 1 f 2

które to rozwiązanie przyjmuje się do dalszych obliczeń.

PRZYKŁAD 4. V przykładzie tym podaje się dwie metody jed- nokrokowe ^ d r u g i e g o rzędu, w których ostatni wzór jest typu (24) oraz takie, że metody , <§2 mają wspólne węzły, co pozwala zredukować liczbę obliczeń wartości funkcji f(x,y).

Kolejno realizuje się wzory s

2 3 h „3

yn+l/4 = yn + Z fn •

3 3 h 2

yn+l/2 = yn + 2 n+1/4 »

yn+1 yn + ^ n + 1 / 2 *

y 3y n + 1

3 -3

Dla wartości y .j i 3^+1 zachodzą związki

(35)

y ( « h ) i y3+, - ^ y 3 (x)h3 y(x+h) i y 3+1 + y 3 (x)h3 .

Została opracowana ^procedura o nazwie diffsysthek (w języku ALGOL 60), która realizuje powyższe wzory. Wyniki obliczeń przedstawiono w § 5 tej pracy. Realizacja metod wyma- ga ^ obliczeń wartości funkcji f (x,y) w jednym etapie całkowa- nia. Krok całkowania jest wyznaczony metodą podstawową, z tym że korzysta się z liczby cń określonej w (29 )• CJ

Konstrukcja metod jednokrokowych ^ dla wyższych rzę- dów jest na ogół trudna. Poniżej podano przykład metody ,

^ 2 trzeciego rzędu.

PRZYKŁAD 5.

3 3 1

yn+2/9 “ + 9

(36)

M.SZYSZKOVICZ y5 i = n+1 Jn o \ n + J h(flł + kr** ,n+1/2 + f** )n+1 '

¥ tym wypadku = Yq^, f2 = ® ’ tzn. rozwiązanie y^+1 ma lo- kalny rząd O(h^). Dla wartości y^+1 i y^+1 nie wykonuje się ekstrapolacji Ricbardsona* Do dalszycłi obliczeń przyjmuje się y5 Jn+1^{i •} O

Metoda trzecia

Jako metody jednokrokowej używa się metody Bobkowa $ rzędu p, której ostatni wzór jest postaci

(30) y (x+h,h) = y £ \ = yP+2+h(A0fP+2+ £ AifP+ ).

1=1 i

V trakcie wyznaczania wartości y(x+h,h/k ) przy użyciu metody *§

realizowany jest wzór postaci

C3’ ) T (“ I-I ) - C U = yn+2łl ±§ ).

a następnie wzory

(32) y(xłfr+1) )= ^ J r+1)h/k =

= ynlr/k + l(Vntr/k + J V AifLr/kw^/k } (r = 1j 2j « « « | k— 1 )•

(37)

Przy odpowiednich węzłach Ci= °» !»•••» <ł ) można po wyzna- czeniu wartości y(x+h,h/k) skonstruować wartość y(x+h,h) bez obliczania funkcji f(x,y). W tym celu wystarczy wykorzystać wartości obliczone przy wyznaczaniu wartości y(x+h,h/k)*

Przyjmuje się następujące oznaczenia:

f k - r / W l + 2 ! i - a (2) k - a i > k + k ~ i » •••» *Z± + k + k “ i- a (k;

Ci = 0, 1; • • • , q )«

Jeżeli w zbiorze węzłów (i = 0, 1,..., q)(r=1, 2, k) występują wszystkie węzły (i = 0, 1,..., q ), to rozwiązanie y(x+h,h) można skonstruować na podstawie wyznaczo- nych wcześniej (tj* podczas obliczania y (x+h,h/lc)) wartości

^n+oć. * ^ przypadku gdy wśród węzłów występują inne niż w {cc.(r)} t dla takich węzłów trzeba obliczyć wartości •

PRZYKŁAD 6. Dana jest metoda Bobkowa rzędu 3

yn + l A = yn + S 3 h h .2

yn+1/2 3 yn + 2 n+1/4 >

3 4

yn+1 *" yn n+1/2 1

(38)

166 M. SZTSZKOVXCZ

Węzły -jot^ j to: węzły oraz z k = 2

to: {°»T5>*£ }i {^ł^ł1 } • A więc w zbiorze węzłów {oL^1^ }, r=1 , 2, występują wszystkie węzły Ot^, zatem rozwiązanie y(x+h,h} można skonstruować bez obliczania wartości funkcji f(x,y), CU

Koszt realizacji metody podstawowej doboru kroku całkowania można obniżyć powyższym sposobem do kl obliczeń wartości funk- cji f(x,y). Korzystając z metody pierwszej obniżania kosztu obliczeń funkcji f(x,y), można zamiast wzoru (3 0) stosować wzór

i podobnie zamiast (31 ) i (32 ) stosować wzory

Użycie powyższych wzorów obniża koszt realizacji metody pod- stawowej doboru kroku całkowania do k(l-1 ) obliczeń wartości funkcji f(x,y).

(39)

Metoda czwarta

Metody jednokrokowe Bobkowa mają często w naturalny sposób charakter metod składanych (embedding method). Tak jest) jeże- li węzeł OL^ = 1• Wtedy otrzymuje się dwie wartości rozwiąza- nia przybliżonego w punkcie x+h

<łi

(33 ) _i y^ . = yp+1 + b ^n+l •'n \ B. fp ^i

oraz

(3^) ' Jn+1 yP+J = yP+1 + b £ A,fP „ .n i n+0^

Parametry {^b± }, C1 = °» Qj) we wzorze (33) spełniają równania, które zapewniają otrzymanie wartości y^+^ • Ola rozwiązań otrzymanych w (33) i w (3*0 zachodzą związki

WYZNACZANIE KROKU CAŁKOWANIA W METODACH 167 JEDNOKROKOWYCH B O B K O W A ___________________________________

y (*ta) S + C2 y (P+,)(x)hP+1

Oznaczając przez H krok całkowania, który daje przybliżone rozwiązanie z zadaną dokładnością oraz przez 8 wektor błędów dopuszczalnych, otrzymuje się nierówność

e * | y (x+H)- y £ ]

Podobnie jak w metodzie podstawowej krok H jest szukany w po—

(40)

Określając

i przyjmując, że zachodzi równość przybliżona

|A'|* |A| J s L kp-1

dla ustalonej wartości k, gdzie A jest wielkością określoną w (10), można podać następujący wzór na liczbę 10:

Przyjmuje się, że const 1. Wartość const ustalana jest do- świadczalnie. Podobnie postępuje się w metodach Rungego-Kutty-

— Fehlberga•

Obliczenia wykonuje się według schematu blokowego zamiesz- czonego w § 2 tej pracy. Zmienione są tylko nieznacznie kry- teria akceptacji rozwiązań przybliżonych » y^+] • ^ak<> roz- wiązanie i wartość do dalszych obliczeń przyjmuje się yj^*]•

PRZYKŁAD 7« Przykładem metody składanej Bobkowa czwartego rzędu jest metoda określona wzorami

(41)

I WYZNACZANIE KROKU CAŁKOWANIA W METODACH

(f5 - 3f3 1 /«, + \ n n+1/3 bf**n+1/2' 1 )

gdzie m jest parametrem i m > 0 . Metody określone w przykła- dzie 5 i 6 są również metodami składanymi. □

5. ILUSTRACJE NUMERYCZNE

W tym miejscu zostaną zilustrowane tylko niektóre z powyżej prezentowanych metod doboru kroku całkowania. Pełne zestawie- nie wyników obliczeń ze stosowania wszystkich podanych w tej pracy metod doboru kroku, zostanie przedstawione w przygotowy- wanym aktualnie raporcie Instytutu Informatyki Uniwersytetu Wrocławskiego. W raporcie tym zamieszczone>również będą teks- ty procedur (w języku ALGOL) realizujących poszczególne meto-

W tej pracy jako zadanie testowe przyjęto zagadnienie po- czątkowe postaci

dy.

y 1 = 1/y2* y i ) = 1 • y -|Cx ) = ®X » y f2 = -l/yt, y 2( 0 ) = 1 , y 2 (x ) = e“X .

(42)

Rozwiązania są wyznaczane kolejno w punktach x = .5, 1.0, 1.5, 10.0 z następującymi dokładnościami eps = ^ “3» 10~9 Wyniki przedstawiono tylko dla x = 1.5 i x = 10.0. Jako wynik podana jest wartość globalnego błędu otrzymanego rozwiązania « (błąd względny) oraz (f J - liczba obliczeń funkcji f(x,y).

V tabeli 1 przedstawiono wyniki realizacji metody podsta- wowej doboru kroku całkowania dla metody Rungego-Kutty czwar- tego rzędu. Obliczenia przeprowadzono z k = 2, 3» 4 i wyniki podano dla x = 10.0.

Tabela 1

eps (1:2) C f J (1:3) i f } (1:4) rfJ

io“3 3^ 1 0-3 -‘‘•7,0-3

10-9 ' ^k T 9 2016 2-°10-9

V przypadku stosowania metody podstawowej doboru kroku całkowania z k = 2, 3» 4 do innych zagadnień testowych, po- twierdza się możliwość stosowania tych podziałów kroku cał- kowania w praktyce.

V tabeli 2 i 3 zestawiono wyniki otrzymane z realizacji metody jednokrokowej drugiego rzędu. W tabeli 2 przedstawio- ne są wyniki otrzymane ż jawnej metody trapezów realizowanej z metodą podstawową doboru kroku całkowania, a w tabeli 3 wy- niki otrzymane z realizacji metody jednokrokowej z przykładu 4 tej pracy.

1.210-2 144 1«910-2 180 -1.8 10-2 2.61q-2

-6.4io~9 1808 -8.41q-9 1700 7-9 10~9 1,010"8 O

(43)

Tabela 2

X 10~3 [ f ] i o " 6 [ f l 10“ 9 [ f ]

1 . 5 - 1 * 11 0 - 4 19 - 7 .2i o- 8 ¹¹⁴ - 1 . 8 , o - 10 1089

6 - 9 i o " 3 - 1 . 7 i o - 7 4 . 2 10- 1°

1 0 . 0 1 . 4i q- 2 203 1 * 5 1 o “ 3 1858 t . 6 , 0 - 8 18443

“ 1 , 8 t o “ 2 - 2 . ° 1 0 - 5 - 2 - 1i o- 8 □

Ta be la ³

X 10~3 [ f ] i o - 6 [ f ] 1 0 ' 9 [ f ]

1 . 5 “ 2 # 7 10” ił 4 - 1 ^ 1 0 - 7 31 5 * 1 1 0 ~ 11 225

2 - 3 10“ 4 1 •3 10” 7 9 . 7 10- 1 1

1 0 . 0 - 2 *3 10‘ 2 54 " 2 *5 10“ 6 4 4 2 -6.2i o -1° 4 2 6 6

2 - ° i o - 2 2 . k i 0 - 6 6.6 1 0 - 1 0 ^D

Podobnie zostały zestawione wyniki w tabeli 4 i 5. V ta- beli 4 przedstawione wyniki obliczeń otrzymano standardową metodą Rungego-Kutty czwartego rzędu z metodą podstawową do- boru kroku całkowania (k = 2), a w tabeli 5 podano wyniki z realizacji metody Bobkowa czwartego rzędu (wzory tej metody są podane w przykładzie 7)» gdzie krok całkowania dobiera się według metody trzeciej doboru kroku całkowania dla metod Bobkowa•

(44)

T a b e l a z

k

1 0 - 3 C f J _{1 0 - 6} i t ] _{1 0 ~ 9} _{C f ]}

1 . 5 9 . 9 , 0 - 7 1 2 6 . 8 i o - 9 4 8 2 . 5 1 0 “ 11 ^{3 2}

1 , 2 1 0 - 6 2 - 2 1 0 " 8 - I .^{6 1 0}- H

1 0 . 0

3 ^ 1 0 - 3 1 5 6 - 2 . 2 i o “ 6 5 2 8

“ 1 * 6 1 0 “ 9 ^{2 0 1 6}

- i,' 7 1 0 - 3 2 . ? i o - 6 2 - ° i o - 9 ⁰

T a b e l a 5

X 1 0 - 3 C f 3 1 0 " 6 C f J 1 0 “ 9 C f j

1 . 5 6 . 3 , 0 - 6 9 1 . 2 i o - 6 2 6

6 ^ 1 0 - 9 6 2

- 3 . 3 ,^{0 - 6} - 7 . 7 i o - 7 - 5 . 7 i o - 9

1 0 . 0 - 1 . 7 ,^{0 - 2} 7 9 - 1 . 7 l 0 - 5 2 2 3 Os . CO ₀ 1 00 8 1 7

2 * 8 , o - 2 2 - 8 1 0 - 5 - 5 . 0 1 0 - 8 a

PRACE CYTOWANE

£ l J J.C h o m i c z, A.O l e j n i c z a k , M.S z y s z k o — w i c z, A method for finding the step size of integration of a system of ordinary differential equations. Zast.Mat.

17 (1983), 645-654.

C2J V.X.K r y 1 o v, V.V.B o b k o v, P.I.M o n a s t y r - n y j, Vyclclitel*nyei metody, t.I, Moskva 1976, t.II, Moskva 1977.

Wyznaczanie kroku całkowania w metodachjednokrokowych Bobkowa*

%)[%f

c±w i%r' = ci ( x+t ) ( i r ' = ••• = «s1(»(^ < > i)(ijp+1

£li.

(

)

j

A

= 1 ■^‘iP^'i = "7+7 ^ =

- 75* i'<i •

- 75* i'<i ^•