Grażyna HOBOT (Lublin)
O pewnej metodzie całkowania równań zwyczajnych z efektywnie różniczkowalną prawą stroną 1. R'ozważmy równanie różniczkowe
(1.1) x' = f (t, x)
z warunkiem początkowym
x(to)=xo.
W pracy tej opisane1 b«idą dwie, analogiczne do metod Rungego-Kutty, metody jedno- krokowe rzCi,du czwartego.
Przy założeniu, że istnieje możliwość praktycznego obliczania dla dowolnie zadanych war·
tości t i x wartości funkcji
(1.3) g t, ( ) _a/(t,x) a/(t,x)/( ) X - at + ax t, X '
przybliżoną wartość rozwiązania problemu ( 1.1) - ( 1.2) w punkcie t k + 1 przedstawimy w po- staci
(1.4)
Funkcja~ jest kombinacją. liniową wartości funkcji/ i g w pewnych punktach pośred
nich.
2. Pierwszą. metodc; jednokrokową rzCidu czwartego określimy w sposób następujący:
(2.1) gdzie
g0 = 0.5h2g (tk, xk), ko = hf (tk, xk),
g1 =0.5h2g(tk +M1h, xk +L1oko +P1ogo), k1 =hf(tk +M1h, xk +R1oko +E1ogo +E11gt)·
1 Pełne dowody przedstawionych tutaj tez zostaną opublikowane osobno.
38 G. Hobo t
Współczynniki ao, a1, bo, bi, L 10, P1o , R 10, E 10 i E 11 zostały wyznaczone w podobny sposób jak w metodzie Rungego-Kutty (patrz Ralston [2], str. 203), a mianowicie z warun- ku, aby metoda była rzędu czwartego i tak dobrane, że rozwiązanie w punkcie pośrednim t k + _M 1 h zgadza się z dokładnym aż do wyrazu z h 2 .
Ządania te doprowadziły do wyznaczenia stałych w zależności od wolnego parametru M 1 .
Są one odpowiednio równe:
(2.3)
2M~ - 2M1+1
ao = 3
2M1 6Mi - 8M1 + 3 bo = ---6,..._M .... i.---'
2M1 -1 a -1 - 2M~ ' b 1 = 3 -4M1 6Mi .
L1o=R1o=M1, Eio =-3-, 2Mi
Natomiast drug'i z metod określają wzory:
(2.4) gdzie
k0 = hf (tk, xk)'
g1 = 0.5h2g (tk + M1 h, xk + L 10k0), (2.5) k1 = hf (tk + M1 h, xk + R10k0 + E11g1 ),
g2 =0.5h2g(tk +M2h, xk +L2oko +L21k 1), k2 =hf(tk +M2h, xk +R20k0 +R21 k1 +E22g2 ).
Postępując podobnie jak w przypadku poprzednim można \Vyznaczyć parametry wystę
pujące we wzorach (2.4) - (2.5). Przeprowadzone obliczenia da.ty
M = 3 -4M1
2 2 (2 - 3Mi)' 6M 1 M2 - 3 (M 1 + M2) + 2
ao = 6M1M2 ' 3M2 -2
L1o=M1,
(2.6) 2
2 3 2
R _ M2 (M2 - M1 + 8M1M2 - l8M1M2 + 6M1M2 + 6M1 - 4M2)
20 2M1 (2M2 -Mi)(2 - 3Mi) '
M2 (M2 -Mi) [4(M2 +Mi)-6M1M2 -1]
R21 =~----~~~---=--~~~~----~~~~~
2M1 (2M2 - Mi)(2 - 3Mi) ' M2 (M2 - M1 - 3MiM2 + 4Mi - 2M1M2) E22 = (2M2 - M1 )(2 - 3Mi)
L21 = M~ ?M ..
- 1
3. Pa~ametr M1 występują.cy we wzorach (2.1) - (2.3) można dobrać w ten sposób, aby uzyskać możliwie małe oszacowanie współczynnika as przy h5 w rozwini<tciu xk+l - x(tk+l) w szereg potęgowy względem h.
Współczynnik ten można przedstawić w postaci (3.1)
8
as= L diei = (d (Mi), e).
i=l
Składowe wektora e (pochodne cząstkowe wyższych rz<tdÓw funkcji f w punkcie t k , x k) są liczbami, natomiast składowe wektora d są wielomianami (stopnia co najwyżej drugiego) ze
wzgl~du na M 1 •
Z (3.1) mamy
(3.2) las I~ ld(M1 )I !el
Długość d (M 1 ) bttdzie najmniejsza dla
(3.3) M1 = 0.6403744628.
Przy takim doborze M 1 uzyskujemy najmniejsze możliwe oszacowanie dla wielkości las I
(co nie oznacza oczywiście, że as osi<Aga minimum). .~
Podstawiając wartość M1 z (3.3) do (2.1)-(2.2) otrzymujemy
(3.4) xk+l = xk + 0.46545275k0 + 0.534547242k1 + 0.0981256692g0 + 0.2172535262g1 , gdzie
k0 = hf(tk, xk), g0 = 0.5h2g (tk, xk),
(3.5) g1 = 0.5h2g (tk + 0.6403744628h, xk + 0.6403744628k0 + 0.4100794514g0 ), k1 = hf (tk + 0.6403744628h, xk + 0.6403744629k0 + 0.2733863006g0 +
+ 0.1366931505gl ).
Przyjmując, że znane są ograniczenia
(3.6) lf1 <M, i+ j ~4.
dostajemy
(3.7) las I~ 0.0694585886L 4 M.
Dla M 1 = 0.5 otrzymujemy jako przypadek szczególny wzory podane przez Zurmiihla l 3 ]
(3 .. 8) gdzie
40
(3.9)
W metodzie tej ' (3 .10)
G. Hobo t
go = 0.5h2g (tk, xk), ko= hf (tk, xk),
g1 = 0.5h2g (tk + 0.5h, xk + 0.5k0 + 0.25go).
4. Metoda (3.4) - (3.5) lepiej reaguje na nagłe zmiany krzywizny rozwiązania w por6w- naniu z klasyczną. metodą Rungcgo-Kutty rzli,du czwartego opisaną wzorami
(4.1) gdzie
(4.2)
ko = hf (tk, xk),
k1 =hf(tk +0.5h, xk +0.5k0), k2 = hf (tk + 0.5h, x,,, + 0.5k1 ), k3 = hf (tk + h, xk + k2).
Natomiast metoda (2.4)-(2.6) daje lepsze wyniki w przypadku, gdy rozwią.zanie ma prze- bieg bardzo spokojny.
Obliczenia testowe były przeprowadzone na maszynie cyfrowej Odra-1013 ze stałym krokiem obliczeń. Otrzymane wyniki por6wnywane były z metodą. opisaną. wzorami (4.1).- (4.2).
W obliczeniach wzi~to pod uwagCi wzory (3.4)-(3.5) oraz (2.4)-(2.6) z tym„ że w tych . h . M l
ostatnie przyJCito 1 = 3.
J. X 1 = X + t + l, x(-1) = 0, h = 0.1;
rozwiązanie dokładne: x(t) = exp (t+ 1) - 2 - t.
tk xk z (3.4)-(3.5) błą.d xk z (2.4).:.-(2.6) błąd -0.60 .918246445610-l -.5310-7 .918248799210-l .1810-6 -0.10 .559602914010+0 -.2010-6 .559603786810 +O .6610-6 0.50 .198168897110+1 -.6010-6 .198169146110+1 .2410-5 1.50 .86824912081o+1 -.2710-5 .868250322410+1 .9310-5 2.00 .160855315310 +2 -.5310-5 .160855552710 +2 .1810-4 tk xk z (4.1)-(4.2) błą.d
-0.60 .918242400810-l -.4610-6 -0.10 .559601412410 +O -.1710-5 0.50 .198168391110+1 -.5110-5 1.50 .868247061610+1 -.2310-4
I 2.00 .160854907510+2 -.4610-4
I (1)1
II. x = -x ctg t ? , x(l). = 1, h = 0.1;
rozwiązanie dokładne: x(t) = sin l(l) sin(-}).
tk xk z (3.4)-(3.5) błąd xk z (2.4)-(2.6) błąd 1.10 .93757824210+0 -.5910-7 .93757888010+0 .4710-8 1.50 . 73486666610+0 -.9110-7 . 73486760810+0 .3210-8 2.60 .44588838210+0 -.5710-7 .44588895 710+0 .4010-9 3.00 .38883604810+0 -.4910-7 .38883654610+0 .4010-9
tk xk z (4.1)-(4.2) błcid 1.10 .93757917410+.o .3410-7 1.50 . 73486808910+0 .5110-7 2.60 .44588927910+0 .3310-7 3.00 .3888-3682810+01 .2910-7
Całkowanie ia pomocą podanych wzorów wymaga oprócz podprogramu obliczającego wartości funkcji/ także podprogramu obliczają.cego wartości funkcji g. Czas wykonywania 'tego ostatniego bywa dłuższy od czasu wykonywania podprogramu obliczającego wartości
funkcji/i w związku z tym czas obliczeń może si€i wydłuiyĆ.Jest to rekompensowane
zwiększoną dokładnością. otrzymanych wyników.
Czasami zdarza się, że funkcja g ma prostszą postać niż funkcja/. Ma to miejsce np.
w przyp~dku równania
1 t3 (4+3t) X
X= l+t +l+t'
dla którego g(t, x) = 12t2 • W takich przypadkach czas obliczeń zmienia si~ w sposób niezna- czny.
Literatura
[ 1] G. H o b o t, Pewne metody rozwiązywania równań różniczkowych zwyczajnych, praca doktorska UMCS Lublin, 1971 (w przygotowaniu do druku).
(2) A.Ra 1 st o n, Wstęp do analizy numerycznej, Warszawa 1971.
(3) R. Z ur m ii h.l, Runge-Kutta Verfahren unter Vervendung hoherer Ableitung, Z. Angcw. Math.
Mech. 32 (1952), str. 153-154.
