Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 1
Rozwiązywanie równań różniczkowych zwyczajnych
Równaniem różniczkowym zwyczajnym n-tego rzędu będziemy nazywać równanie postaci:
dy d2y d3y dny
F x,y, , , , , 0
2 3 n
dt
dx dx dx
Ogólnie równanie takie może być spełnione przez wiele funkcji y (x )  aby możliwe było jednoznaczne wyznaczenie
rozwiązania konieczne jest podanie dodatkowych informacji  wartości funkcji y (x) i ewentualnie jej pochodnych dla
jakiejś wartości x. Koniecznych będzie n takich wartości. Jeśli podano je wszystkie dla tego samego punktu x to mamy
0
zagadnienie początkowe, jeśli wchodzi w grę więcej niż jedna wartość x to zagadnienie takie nazywamy brzegowym.
Zazwyczaj w przypadku równań różniczkowych zwyczajnych wyższych rzędów możliwe jest dokonanie podstawienia
n  1 nowych zmiennych (za wyższe pochodne) i wówczas równanie wyższego rzędu zostaje zamienione na układ n równań
rzędu pierwszego. Dlatego też omawia się zwykle tylko metody rozwiązywania równań rzędu pierwszego.
Rozwiązania numeryczne, które będziemy otrzymywać nie będą oczywiście w postaci funkcji analitycznych  będą
to wartości zbliżone do wartości rozwiązania obliczone dla pewnej skończonej liczby punktów x . Odległości pomiędzy
i
dwoma sąsiednimi punktami będziemy oznaczać przez h (na razie zakładamy, że są one stałe, możliwość zmiany kroku
będzie omawiana oddzielnie).
Zagadnienia poczÄ…tkowe
Rozwinięcie w szereg Taylora
Jedną z metod rozwiązywania równań różniczkowych zwyczajnych postaci dy dx f x,y jest wyrażenie
rozwiÄ…zania y (x) w postaci szeregu Taylora w otoczeniu punktu x0:
2 3
h h
y x0 h y x0 h f x0,y x0 f x0,y x0 f x0,y x0
2! 3!
Jeśli jako warunek początkowy podano y (x0) to f (x0,y (x0)) otrzymamy bezpośrednio z równania różniczkowego, a
pozostałe pochodne można wyliczyć. Pożądaną ilość wyrazów należy wstawić do badanego równania i scałkować. Niestety
w ogólnym przypadku obliczenia kolejnych pochodnych stają się zbyt skomplikowane.
Metoda Eulera
Jest to najprostsza (w swej klasycznej postaci) z metod numerycznych. Odpowiada ona rozwiÄ…zaniu za pomocÄ…
rozwinięcia w szereg Taylora, w którym zachowano tylko pierwsze pochodne. Rozwiązywane równanie różniczkowe
przekształcimy do postaci:
dy
f x,y
dx
i założymy, że podano nam początkową wartość funkcji y0(x0).
Działanie metod Eulera (podstawowej i jej obu modyfikacji) będzie ilustrowane na dwóch przykładowych równaniach
różniczkowych:
dy 2x dy
y i x y
dx y dx
Równania te będą rozwiąywane wszystkimi trzema metodami w przedziale <0,1>, dla warunku początkowego y(0) = 1
i z trzena różnymi krokami h wynoszącymi 0.2, 0.1 i 0.05.
Klasyczna metoda Eulera
Możemy zastosować następujący algorytm:
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 2
y1 y x0 h f x0,y x0
yi 1 yi h f xi ,yi i >0
za pomocą którego wychodząc od znanych wartości początkowych x0,y0 można obliczyć kolejno wartości szukanej funkcji
we wszystkich dalszych punktach x . Można udowodnić, że błąd (zdefiniowany jako różnica rozwiązania analitycznego
i
i obliczonych wartości przybliżonych) dąży do zera dla h dążącego do zera. Metoda Eulera jest metodą zbieżną o
dokładności rzędu O(h) (chodzi tu o błąd globalny  tzn. po przejściu wszystkich kroków  błąd lokalny 
odpowiadający przejściu od jednej wartości do sąsiedniej  jest rzędu O(h 2)).
Zmodyfikowane metody Eulera
Pierwsza modyfikacja metody Eulera polega na obliczeniu najpierw wartości pośrednich:
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 3
h
xi 1 xi
2
2
h
yi 1 yi f xi ,yi
2
2
fi 1 f x ,yi 1
1
i
2 2
2
a następnie wykorzystaniu ich do obliczenia kolejnej wartości y :
i+1
yi 1 yi h fi 1
2
Druga modyfikacja metody Eulera (zwana też metodą Heuna) polega na obliczeniu najpierw wartości  wstępnie
przybliżonych :
yi 1 yi h f xi ,yi
fi 1 f xi 1,yi 1
po czym oblicza się wartości y :
i+1
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 4
fi fi 1
yi 1 yi h
2
W obu modyfikacjach odrzucona reszta rozwinięcia w szereg Taylora jest rzędu O(h3) (jest to błąd lokalny  tzn.
odpowiadający przejściu pomiędzy dwiema kolejnymi wartościami przybliżenia).
W praktyce przyjmuje się, że błąd uzyskanego rozwiązania jest równy różnicy pomiędzy rozważanym rozwiązaniem,
a innym rozwiązaniem tego samego zagadnienia otrzymanym dla dwa razy większego kroku h.
Błąd rozwiązania można też szacować na podstawie tzw. kryterium Rungego  wynosi on różnicy rozważanego
rozwiązania i rozwiązania otrzymanego dla kroku dwukrotnie większego.
Porównanie metod Eulera
Poniżej znajdują się wykresy przedstawiające rozwiązania oraz wartości błędów rozwiązań otrzymanych
poszczególnymi metodami dla różnych długości kroków.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 5
Jak widać nie można jednoznacznie stwierdzić, która z modyfikacji metody Eulera jest lepsza  zależy to od
konkretnego rozwiÄ…zywanego zagadnienia.
Metoda Eulera z iteracjÄ…
Metodę Eulera można udokładnić stosując dla każdej obliczonej wartości y dodatkowo proces iteracyjny. Traktując
i
wartość obliczoną klasyczną metodą Eulera jako zerowe przybliżenie:
yi 0 yi h f xi,yi
1
zastosujemy następujący proces iteracyjny:
h
yi k yi f xi,yi f xi 1,yi k 1
1 1
2
który będzie powtarzany, aż dwa kolejne wyniki będą zgodne z zadaną ilością cyfr znaczących, wtedy ostatni z nich
przyjmiemy za y .
i+1
Zazwyczaj proces taki jest szybkozbieżny (tzn. zadowalający wynik powiniśmy otrzymać już po kilku iteracjach) 
jeśli po 3-4 iteracjach nie uzyskujemy zbieżności należy zmniejszyć krok h.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 6
Rozwiązywanie równań ruchu
Bardzo częstym przypadkiem jest rozwiązywanie równań ruchu (dla jednego lub bardzo wielu ciał)  są to oczywiście
równania drugiego rzędu. Każde takie równanie (drugiego rzędu) można zapisać jako układ dwóch równań pierwszego
rzędu. Poniższe równania zostały zapisane wektorowo  z punktu widzenia rozwiązania numerycznego liczba
występujących w problemie wymiarów sprowadza się do wymiaru tablic wykorzystywanych do przechowywania
rozwiÄ…zania.
Podstawową postać równania ruchu:
d2R F
a
2
m
dt
Przepiszemy w nieco innej postaci:
dv
a r,v
dt
dr
v
dt
Typowy algorytm Eulera będzie mieć wówczas taką postać:
vi 1 vi Äai
ri 1 ri Ävi
Algorytm Eulera-Cromera
PewnÄ… modyfikacjÄ… algorytmu Eulera jest algorytm nazywany algorytmem Eulera-Cromera. Modyfikacja polega na tym,
że do obliczenia kolejnej wartości położenia wykorzystuje się obliczoną przed chwilą wartość prędkości:
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 7
vi 1 vi Äai
ri 1 ri Ävi 1
Jest to bardzo niewielka zmiana (co do zasady identyczna z modyfikacją Seidla, którą poznaliśmy w przypadku układów
równań), jednak w większości przypadków algorytm Eulera-Cromera zdecydowanie góruje nad typowym algorytmem
Eulera.
Algorytm Verlet a
Załóżmy, że mamy do rozwiązania równanie ruchu (albo równania  dla wielu ciał). Ponieważ jest to równanie 2-go
stopnia rozwiązujemy je jako układ dwóch równań rzędu pierwszego. Jeśli interesuje nas przede wszystkim trajektoria
ciała, to można byłoby twierdzić, że rozwiązanie pierwszego z równań  dające prędkość  jest dla nas stratą czasu.
W sytuacji takiej zastosowanie znajduje algorytm Verlet a.
Jeszcze raz równanie ruchu:
2
d R F
a
2
m
dt
Drugą pochodną zapiszemy w postaci różnicowej:
d2R 1
R 2R R
n 1 n n 1
2
dt Ä2
gdzie n oznacza numer kolejny kroku czasowego, a Ä jest jego dÅ‚ugoÅ›ciÄ….
Proste przekształcenie doprowadza nas do:
R 2R R Ä2a
n 1 n n 1 n
Jak widać z wzoru na Rn+1 do  uruchomienia tego algorytmu potrzebne są dwie kolejne wartości początkowe położeń.
W problemach fizycznych z jakimi możemy się spotkać, zazwyczaj będziemy znać tylko jedną wartość położenia  w
chwili początkowej i odpowiadającą jej wartość prędkości. Dla  uruchomienia algorytmu konieczne będzie
 wygenerowanie drugiej potrzebnej wartości  najprościej można to zrobić zakładając stałość przyspieszenia w
pierwszym odcinku czasu Ä. W takim przypadku wartość R1 bÄ™dzie taka jak dla ruchu jednostajnie przyspieszonego:
Ä2
R1 R0 ÄV0 a0
2
Gdyby potrzebna była również prędkość cząstki to ją również możemy zapisać w postaci różnicowej (skorzystamy
z wzoru centralnego):
dR 1
V R R
n 1 n 1
dt 2Ä
skÄ…d otrzymamy:
R R
n 1 n 1
V
n
2Ä
Jak widać wzór ten może w pewnych sytuacjach być  nieporęczny  obliczana prędkość jest  opóz
niona w
stosunku do położenia o jeden krok czasowy. Można temu zaradzić obliczając ją wprost z działającej siły:
V V Äa
n 1 n n
znacznie dokładniejsze wyniki uzyskamy wykorzystując do obliczeń zamiast an średnią wyznaczoną z an i an+1 (an+1
będziemy znać bo  przed chwilą obliczyliśmy położenie Rn+1).
Poniżej przedstawiono rozwiązania prostego zagadnienia: harmonicznego ruchu wahadła bez oporów, tłumienia i
wymuszenia. Dragjąca masa jest jednostkowa, współczynnik sprężystości też jest jednostkowy.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 8
Pierwsze trzy rysunki przedstawiają wykres położenia drgającej masy w funkcji czasu. W przypadku podstawowej
metody Eulera już  na oko widać, że wyniki są złe  wiemy, że energia, a co za tym idzie i amplituda, opisywanego
układu powinny być stałe. Jednak wybór którejś z pozostałych dwóch metod może być trudny. Aby lepiej je ocenić
wyznaczymy całkowitą energię układu i zbadamy jej zmiany z czasem. Przy rzyjętych parametrach (jednostkowa masa
i współczynik sprężystości) całkowita energia wynikająca z warunków początkowych również jest jednostkowa. Fakt, że
energia ta ewidentnie nie jest stała podczas trwania ruchu, nie musi od razu dyskwalifikować metody  ze względu na
różnego rodzaju błędy nie możemy oczekiwać rzeczywistej stałości całkowitej energii  decydującym kryterium będzie
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 9
to czy obliczona energia całkowita waha się w ograniczonym zakresie czy też amplituda jej wahań wzrasta. Oczywiście
celem wyboru metody należałoby sporządzić te wykresy dla dużo dłuższego czasu  poniżej czas przedłużono
dziesięciokrotnie:
Metody Rungego-Kutty
Metodę Eulera otrzymujemy dzięki rozwinięciu szukanej funkcji w szereg Taylora i odrzuceniu wyrazów
zawierających pochodne rzędu wyższego niż pierwszy  błąd lokalny jest wówczas rzędu O(h2). Zachowanie wyższych
pochodnych jest możliwe tylko dla bardzo szczególnych postaci równań, w ogólnym przypadku będą one zbyt
skomplikowane. Możliwe jest jednak stworzenie metody jednokrokowej (tzn. wymagającej znajomości tylko jednego
wcześniejszego punktu) o dokładności lokalnej wyższego rzędu. Są to metody Rungego-Kutty.
Ogólnie algorytm metod Rungego-Kutty jest następujący:
yi 1 yi h Ć xi,yi,h
funkcja Ć(x,y,h) jest odpowiednio dobraną aproksymacją funkcji f (x,y ) na przedziale (x ,x ).
i i+1
Metoda Rungego-Kutty rzędu drugiego
Załóżmy, że funkcję Ć(x,y,h ) zdefiniujemy jako średnią z dwóch przybliżeń pochodnej f (x,y ), oznaczymy je przez
k1 i k2: Ć(x,y,h) = ak1 + bk2. Wówczas algorytm Rungego-Kutty (drugiego rzędu) będzie mieć postać:
yi 1 yi h ak1 bk2
Załóżmy dalej, że:
k1 f xi,yi
k2 f xi ph,yi qhf xi,yi f xi ph,yi qhk1
gdzie p i q są stałymi, które wyznaczymy póz
niej.
Rozwińmy k2 w szereg Taylora opuszczając wszystkie wyrazy zawierające potęgi h wyższe niż pierwsza:
k2 f xi ph,yi qhf xi,yi
2
f xi,yi phfx xi,yi qhf xi,yi fy xi,yi O h
Szereg Taylora dla funkcji dwóch zmiennych:
2 2
r s
f (x r,y s) f (x,y) rfx sfy fxx rsfxy fyy
2 2
Podstawiając definicję k1 i rozwinięcie k2 do algorytmu otrzymamy:
yi 1 yi h af xi,yi bf xi,yi
2 3
h bpfx xi,yi bqf xi,yi fy xi,yi O h
Następnie rozwiniemy w szereg Taylora funkcję będącą rozwiązaniem dokładnym:
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 10
2 3
h h
y xi h y xi 1 y xi hf xi ,y xi f xi ,y xi f ¾,y ¾
2 3!
gdzie ¾ jest liczbÄ… z przedziaÅ‚u (xi,xi+1). Pochodna f (xi,y(xi)):
f xi ,y xi fx xi ,y xi fy xi ,y xi f xi ,y xi
Jeśli teraz zażądamy, aby współczynniki przy odpowiednich potęgach h były równe to otrzymamy:
y xi yi
f xi ,y xi a b f xi ,yi
1
fx xi ,y xi fy xi ,y xi f xi ,y xi bpfx xi ,yi bqfy xi ,yi f xi ,yi
2
Jeśli założymy yi = y (xi) i zażądamy, aby dla wszystkich funkcji współczynniki przy h2 były równe otrzymamy
a + b = 1, bp = ½ oraz bq = ½, a zatem:
a 1 b
1
p q
2b
Ponieważ mamy tylko trzy równania, a cztery niewiadome, możemy jedną z nich wybrać dowolnie, najczęściej
przyjmuje siÄ™ b = ½ albo b = 1.
W przypadku b = ½ otrzymamy a = ½, p = 1 i q = 1. Algorytm bÄ™dzie mieć wówczas postać:
h
yi 1 yi f xi ,yi f xi h,yi hf xi ,yi
2
co można zapisać również w postaci:
yi 1 yi hf xi ,yi
h
yi 1 yi f xi ,yi f xi 1,yi 1
2
Jest to metoda identyczna z podaną wcześniej drugą modyfikacją metody Eulera. Metody tego typu nazywane są
metodami predyktor-korektor, ponieważ najpierw za pomocą metody dającej gorsze przybliżenie przewidujemy wartość
y (jest to wartość jaką yi+1 miałoby gdyby nachylenie (pochodna) wynosiło k1), a następnie uzyskujemy skorygowaną
i+1
wartość yi+1 (jako średnią ważoną z wartości f na obu końcach przedziału (xi,xi+1). Metoda ta jest znana również pod nazwą
metody Heuna. Ilustruje ją lewy rysunek powyżej.
Dla b = 1 otrzymamy a = 0, p = ½ i q = ½ i algorytm przyjmie postać:
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 11
h
yi 1 yi f xi ,yi
2
2
h
yi 1 yi hf xi ,yi 1
2
2
Ten algorytm z kolei odpowiada przedstawionej już pierwszej modyfikacji metody Eulera. Najpierw za pomocą
normalnej metody Eulera otrzymujemy wartość y w połowie przedziału (xi,xi+1), a następnie znajdujemy wartość f w tym
punkcie i wykorzystujemy ją do obliczenia y na końcu przedziału. Ilustruje to prawy rysunek powyżej.
Metoda Rungego-Kutty rzędu trzeciego
Metody Rungego-Kutty wyższych rzędów wyprowadza się analogicznie. Dla metody rzędu trzeciego funkcja Ć ma
postać: Ć = ak1 + bk2 + ck3, gdzie k1, k2 i k3 są przybliżeniami wartości pochodnej w różnych punktach wewnątrz przedziału
(xi,xi+1):
k1 f xi ,yi
k2 f xi ph,yi phk1
k3 f xi rh,yi shk2 r s hk1
Algorytm jest więc dany wzorem:
yi 1 yi h ak1 bk2 ck3
Aby wyznaczyć stałe a, b, c, p, r i s rozwijamy k2 i k3 oraz funkcję y w szeregi Taylora wokół punktu (xi,yi) odrzucając
wyrazy zawierające potęgi h wyższe niż trzecia i przyrównujemy odpowiednie wyrazy. Otrzymamy układ równań:
a b c 1
1
bp cr
2
1
2 2
bp cr
3
1
cps
6
Dwie z tych stałych można przyjąć dowolnie. Kutta dokonał wyboru, przy którym algorytm dla metody trzeciego rzędu
jest następujący:
h
yi 1 yi k1 4k2 k3
6
k1 f xi ,yi
1 1
k2 f xi h,yi hk1
2 2
k3 f xi h,yi 2hk2 hk1
W powyższych wzorach przyjÄ™to a = c = 1/6, b = 4/6, p = ½, r = 1 i s = 2.
Poniżej przedstawiam błędy metody Runge-Kutty trzeciego rzędu zastosowanej do identycznych jak poprzednio dwóch
równań przykładowych:
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 12
Metoda Rungego-Kutty rzędu czwartego
Algorytm metod rzędu czwartego jest oczywiście następujący:
yi 1 yi h ak1 bk2 ck3 dk4
Używanych jest kilka algorytmów, autorem dwu poniższych jest Kutta:
h
yi 1 yi k1 2k2 2k3 k4
6
k1 f xi ,yi
1 1
k2 f xi h,yi hk1
2 2
1 1
k3 f xi h,yi hk2
2 2
k4 f xi h,yi hk3
albo
h
yi 1 yi k1 3k2 3k3 k4
8
k1 f xi ,yi
1 1
k2 f xi h,yi hk1
3 3
2 1
k3 f xi h,yi hk1 hk2
3 3
k4 f xi h,yi hk1 hk2 hk3
Jednak najczęściej wykorzystywaną metodą czwartego rzędu jest metoda Gilla:
h
1 1
yi 1 yi k1 2 1 k2 2 1 k3 k4
2 2
6
k1 f xi ,yi
1 1
k2 f xi h,yi hk1
2 2
1 1 1 1
k3 f xi h,yi hk1 1 hk2
2 2
2 2
1 1
k4 f xi , h,yi hk2 1 hk3
2 2
Poniżej przedstawiam błędy metody Runge-Kutty czwartego rzędu zastosowanej do identycznych jak poprzednio
dwóch równań:
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 13
Ogólny schemat metod Runge-Kutta
Ogólny wzór podany przez Rungego:
k1 f x,y
k2 f x c2h, y ha2,1k1
s 1
k f x c h, y h a kj
s s s,j
j 1
s
Åš x,y,h biki
i 1
yk 1 yk h Åš xk,yk,h
Odnosi się on w gruncie rzeczy do wszystkich metod Runge-Kutta. Indeks s odpowiada tu rzędowi metody. W
literaturze można spotkać taki zapis zmiennych występujących w tym schmacie:
0
c2 a2,1
c3 a3,1 a3,2
c a a a
s s,1 s,2 s,s 1
b1 b2 b b
s 1 s
Ogólnie rzecz biorąc podczas  własnoręcznego tworzenia jakiegoś schematu tego typu potrzebne będą rozwinięcia
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 14
Taylora funkcji kilku zmiennych. Ogólny wzór na takie rozwinięcie funkcji f w pobliżu punktu a jest następujący:
n
f
f x1,x2, ,x f a1,a2, ,a xi ai
n n
xi a1,a2, ,an
i 1
n n
2
1 f
xi ai xj aj R x1,x2, ,x
m n
2! xi xj a1,a2, ,a
i 1 j 1
n
Błędy metod Rungego-Kutty
Ogólnie błąd w metodach Rungego-Kutty można określić następująco. Załóżmy, że błąd zdefiniujemy jako różnicę
pomiędzy wartością rozwiązania analitycznego y (gdyby było znane), a wartością otrzymaną przy kroku h1 lub h2 w
wyniku scałkowania równania na przedziale (x ,x ):
i i+1
xi xi
m
´1 y y1 Å‚m h1 1 1
h1
xi xi
m
´2 y y2 Å‚m h2 1 1
h2
Dzieląc oba powyższe równania stronami otrzymamy:
h1 m
y1 y2 h2
y
h1 m
1
h2
co dla h2 = h1/2 przechodzi w:
y1 2my2
y
1 2m
Podstawiając tę wartość do pierwszego lub drugiego wyrażenia na błąd otrzymamy:
y2 y1
´1 2m
2m 1
y2 y1
´2
2m 1
gdzie m jest rzędem metody, a y1 i y2 są wartościami otrzymanymi odpowiednio dla kroków h1 i h2 = h1/2. Oczywiście taki
sposób obliczania błędu wymaga praktycznie potrojenia ilości obliczeń, nie nadaje się więc do praktycznego zastosowania.
W przypadku metod nieparzystego rzędu istnieje możliwość przeprowadzenia po każdym cyklu cyklu odwrotnego (tzn.
ponownego obliczenia wartości y na podstawie właśnie obliczonej wartości y ). Po takim obliczeniu jako błąd można
i
i+1
przyjąć połowę różnicy między obliczonymi w obu kierunkach wartościami y .
i
Wybór długości kroku
Collatz podał praktyczne jakościowe kryterium doboru kroku całkowania dla pierwszego z podanych algorytmów
Kutty czwartego rzędu: po każdym kroku należy obliczyć wartość wyrażenia: k3  k2 / k2  k1 . Jeżeli przekracza ona
kilka setnych to długość kroku należy zmniejszyć.
Często stosowana jest metoda polegająca na obliczaniu (co jakąś liczbę kroków) wartości przybliżonej błędu:
´ = y1  y2 / 15. JeÅ›li wartość ´ jest mniejsza od pożądanej dokÅ‚adnoÅ›ci obliczenia prowadzone sÄ… dalej. JeÅ›li jest
mniejsza niż 1/50 założonej dokładności krok powiększa się dwukrotnie. Jeśli jest większa  obliczenia powtarza się (od
miejsca gdzie ostatnio sprawdzano ´) z krokiem dwukrotnie mniejszym.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 15
Najczęściej spotykaną metodą obliczeń ze zmiennym krokiem jest metoda Runge-Kutta-Fehlberg (o której dalej).
W przypadku programu mającego możność zmiany długości kroku musimy oczywiście podać dopuszczalny zakres
tych zmian, a więc podać jego maksymalną i minimalną dopuszczalną długość. Oczywiście wektor wynikowy takiego
programu (zawierający obliczone wartości yi) ma nieprzewidywalną długość. Co więcej, konieczne będzie zapisanie w
pamięci również wektora  chwil czasu , którym odpowiadają wyznaczone wartości funkcji  zatem cały program może
okazać się bardziej  pamięciożerny . Z drugiej jednak strony, jeśli charakter funkcji nie wymusi skracania kroku, a
umożliwi jego wydłużanie, możemy zaoszczędzić i na pamięci i na czasie (a przynajmniej na cyklach obliczeń). Jeśli okaże
się, że uzyskane rozwiązania nie wypadają w punktach, na których nam zależało będziemy musieli przeprowadzić
interpolacjÄ™.
Stabilność metod jednokrokowych
Przykład: Rozważmy znowu równanie:
dy
f x,y x y
dx
którego rozwiązanie analityczne ma postać:
x x0
y x x 1 1 x0 y x0 e
Przy warunku początkowym y (0) =  1 jego rozwiązaniem analitycznym jest: y (x ) =  x  1. Jak widać przy takim
warunku początkowym składnik eksponencjalny w rozwiązaniu analitycznym znika, ale nawet minimalna zmiana warunku
początkowego, np. na  0,99999 spowoduje w końcu znaczną zmianę wartości y dla dużych x. Tymczasem rozwiązywanie
równań różniczkowych metodami jednokrokowymi sprowadza się w każdym kroku do rozwiązywania kolejnego
zagadnienia poczÄ…tkowego  wynik uzyskany w kroku i staje siÄ™ warunkiem poczÄ…tkowym dla obliczenia wyniku i +1.
Tego rodzaju niestabilność jest cechą konkretnego równania i jego warunku początkowego, a nie zależy od tego jaki
konkretny algorytm jednokrokowy zostanie wykorzystany.
Metoda Runge-Kuta-Fehlberg
Metoda ta jest znana również jako  RKF45" (funkcja  ODE45" Matlaba wykorzystuje tę właśnie metodę). Jest to
metoda o zmiennym kroku, jednak w celu doboru długości kroku nie stosuje się w niej porównywania wyników
otrzymanych przy dwóch różnych długościach kroków  zamiast tego porównuje wyniki trzymane dla dwóch różnych
rzędów metody obliczeniowej (czwartego i piątego).
Algorytm.
Kolejne współczynniki:
k1 hf ti,yi
1 1
k2 hf ti h,yi k1
4 4
3 3 9
k3 hf ti h,yi k1 k2
8 32 32
12 1932 7200 7296
k4 hf ti h,yi k1 k2 k3
13 2197 2197 2197
439 3680 845
k5 hf ti h,yi k1 8k2 k3 k4
216 513 4104
1 8 3544 1859 11
k6 hf ti h,yi k1 2k2 k3 k4 k5
2 27 2565 4104 40
Kolejne przybliżenie (obliczenia rzędu 4):
25 1408 2197 1
yi 1 yi k1 k3 k4 k5
216 2565 4104 5
Lepsze rozwiązanie (obliczenia rzędu 5):
16 6656 28561 9 2
zi 1 yi k1 k3 k4 k5 k6
135 12825 56430 50 55
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 16
Następny krok wyznacza się badając następującą wiekość:
1 1
Tol h Tol h
4 4
s 0.84
2 zi 1 yi 1 zi 1 yi 1
Długość kroku możemy regulować na kilka sposobów:
" metodą  konserwatywną : Jeśli s < 3/4  krok zmniejszamy o połowę, jeśli s > 3/2  krok zwiększamy dwukrotnie
 jak łatwo zauważyć przy zastosowaniu tej metody przez większość czasu krok pozostanie stały
" metodą  dynamiczną : Sprowadzamy wartość s do wybranego zakresu (np. <0.1,4.0>  tzn. jeśli s < 0.1 to
podstawiamy s = 0.1 i jeśli s > 4 to podstawiamy s = 4), a następnie mnożymy dotychczasowy krok przez s (postępujemy
tak w każdym kroku).
" metodą  mieszaną : Postępujemy tak jak w metodzie  konserwatywnej , ale zmianę kroku uzyskujemy mnożąc jego
aktualną długość przez obliczoną wartość s.
Oczywiście niezależnie od zastosowanej metody w programie powinny być narzucone granice dopuszczalnych zmian
długości kroku  tzn. maksymalna i minimalna nieprzekraczalna długość kroku.
Podane wartości granic przedziałów wartości współczynnika s wynikają z praktyki, możemy je dobrać do konkretnego
rozwiÄ…zywanego zagadnienia
Istnieje pogląd, że jako wynik procedury należy, zamiast wektora y, przedstawiać obliczany równolegle wektor z.
Poniższe wykresy przedstawiają porównanie metody Runge-Kutta rzędu 4 i metody Runge-Kutta-Fehlebrg
(wykorzystującą przedstawioną wyżej  konserwatywną metodę zmiany kroku) zastosowanych do już wcześniej podanych
dwóch równań przykładowych. Początkowy krok zastosowany w metodzie RKF był taki sam jak krok używany przez cały
czas obliczeń w metodzie RK4. Warto zauważyć, że zastosowanie w metodzie RKF zbyt dużej tolerancji może prowadzić
do znacznych błędów.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 17
Metoda Runge-Kutta-Nyström
Jak wcześniej powiedziano liniowe równania wyższych rzędów możemy zastępować układami równań niższych
rzędów i rozwiązywać je kolejno. Jednak w przypadku równań rzędu drugiego opracowano metodę czwartego rzędu,
będącą adaptacją metody Runge-Kutta. Metoda ta stosuje się do równań postaci:
y f x,y,y
Zakładamy, że znane są wartości funkcji i jej pochodnej w punkcie początkowym. Dla równoodległych punktów xi (krok
równy h) algorytm jest następujący:
1
yn 1 yn h yn k1 k2 k3
3
1
k1 hf x ,yn,yn
n
2
1 1 1 1
k2 hf x h , yn h yn k1 , yn k1
n
2 2 2 2
1 1 1 1
k3 hf x h , yn h yn k1 , yn k2
n
2 2 2 2
1
k4 hf x h , yn h yn k3 , yn 2k3
n
2
1
yn 1 yn k1 2k2 2k3 k4
3
W metodzie tej można stosować opisany już wcześniej dla metody Runge-Kutta sposób zmieniania długości kroku
h w trakcie obliczeń, z tym, że oprócz różnicy wartości y1 i y2 należy kontrolować również różnice pomiędzy
odpowiednimi pochodnymi (i tak przecież obliczanymi); wydłużenie kroku jest możliwe gdy zarówno wartości funkcji
jak i pochodnych spełniają warunek.
Zastosowanie metod jednokrokowych do układów równań
Jak już wcześniej wspominałem, podręczniki nie omawiają sposobów rozwiązywania równań różniczkowych rzędów
wyższych niż pierwszy, ponieważ zawsze można je sprowadzić do układu odpowiedniej liczby równań rzędu pierwszego.
Rozważmy równanie ruchu (czyli równanie rzędu drugiego):
2
r
a
2
t
podstawiając nową zmienną (prędkość) otrzymamy układ dwóch równań (wektorowych) rzędu pierwszego (albo sześciu
równań skalarnych):
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 18
x
v
x
t
y
v
y
t
z
r
v
z
v
t
t
v
x
v
a
a x
t
t
v
y
a
y
t
v
z
a
z
t
W dalszym ciągu będę omawiać skalarne równanie drugiego rzędu, które po przekształceniu do układu dwóch równań
rzędu pierwszego zapiszemy wektorowo:
dr
v
2
dt r v
d r d
a
2
dv dt v a
dt
a
dt
Gdybyśmy rozważali ruch w trzech wymiarach we wzorach powyżej mielibyśmy wektory nie dwu- ale sześcioelementowe.
Aby zastosować metodę Runge-Kutty trzeciego rzędu właściwe dla niej wzory zapiszemy w postaci wektorowej:
ri vi
k1 f ti, yi f ti,
vi a(ri)
1 1
ri hvi vi ha(ri)
ri vi
2 2
1 1 1 1 1
k2 f ti h, yi h k1 f ti h, h f ti h,
2 2 2 2 2 1 1
vi a(ri)
vi ha(ri) a(ri hvi)
2 2
1
ri vi 2ha(ri) vi
k3 f ti h, yi 2h k2 h k1 f ti h, 2h h
vi a(ri 1hvi a(ri)
2
2 1
ri 2hvi h a(ri) hvi
vi 2ha(ri hvi) ha(ri)
2
f ti h,
1
2
vi 2ha(ri hvi) ha(ri)
a(ri 2hvi h a(ri) hvi
2
I ostateczne rozwiÄ…zanie:
h
yi 1 yi k1 4 k2 k3
6
Warto zwrócić uwagę, że do zapisania tych wzorów w Scilabie wyniki przekształceń widoczne po prawej stronie
przedstawionych powyżej wzorów nie będą potrzebne  funkcję f zdefiniujemy jako przyjmującą dwa argumenty, z
których pierwszym będzie skalar (czas), a drugim wektor zawierający wszystkie niewiadome.
Analogicznie będziemy postępować w przypadku każdej innej metody.
Metody wielopunktowe
Wszystkie metody przedstawione powyżej są nazywane jednopunktowymi  a to dlatego, że do prowadzenia obliczeń
wystarczy znajomość tylko jednej wartości szukanej funkcji. Podanie warunku początkowego jest wystarczające. W
przeciwieństwie do nich metody wielopunktowe wymagają znajomości wartości funkcji w kilku punktach początkowych.
Zazwyczaj jednak warunek początkowy zadany jest w jednym punkcie, wówczas możliwe jest rozpoczęcie obliczeń którąś
z przedstawionych powyżej metod jednopunktowych, obliczenie niezbędnej liczby wartości początkowych, a następnie
kontynuowanie obliczeń metodą wielopunktową.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 19
Metoda Adamsa
Załóżmy, że dla zagadnienia y = f(x,y) z warunkiem początkowym y(x0) = y0 znamy trzy kolejne wartości
y1 = y(x1) = y(x0+h), y2 = y(x2) = y(x0+2h), y3 = y(x3) = y(x0+3h)  obliczone w inny sposób.
Rozwiązanie znajdujemy w dwóch etapach: w pierwszym dokonujemy ekstrapolacji znanych punktów i znajdujemy
przybliżoną wartość punktu kolejnego:
h
yi 1 yi 55fi 59fi 1 37fi 2 9fi 3
24
A następnie, wykorzystując obliczoną wartość yi* do obliczenia fi* w etapie drugim dokonujemy jej udokładnienia:
+1 +1
h
yi 1 yi 9fi 1 19fi 5fi 1 fi 2
24
Poniżej przedstawiam błędy metody Adamsa zastosowanej do identycznych jak poprzednio dwóch równań, do
obliczenia punktów początkowych wykorzystano metodę Runge-Kutty drugiego, trzeciego i czwartego rzędu:
Metoda Milne'a
Metoda ta działa analogicznie jak przedstawiona powyżej metoda Adamsa  obliczenia przeprowadzamy również
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 20
w dwóch etapach:
4h
yi 1 yi 3 2fi 2 fi 1 2fi
3
h
yi 1 yi 1 fi 1 4fi fi 1
3
Jest to metoda 4-go rzędu. Przyjmuje się, że błąd bezwzględny wartości udokładnionej (obliczonej w drugim etapie)
jest równy w przybliżeniu 1/29 różnicy pomiędzy wartością uzyskaną w pierwszym etapie a wartością z etapu drugiego
 możliwe jest więc wygodne kontrolowanie dokładności podczas obliczeń.
Poniżej przedstawiam błędy metody Milne a czwartego rzędu zastosowanej do identycznych jak poprzednio dwóch
równań, do obliczenia punktów początkowych wykorzystano metodę Runge-Kutty drugiego, trzeciego i czwartego rzędu:
Istnieje także metoda Milne'a rzędu 6-go:
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 21
3h
yi 1 yi 5 11fi 14fi 1 26fi 2 14fi 3 11fi 4
10
2h
yi 1 yi 3 7fi 1 32fi 12fi 1 32fi 2 7fi 3
45
Poniżej przedstawiam błędy metody Milne a szóstego rzędu zastosowanej do identycznych jak poprzednio dwóch
równań, do obliczenia punktów początkowych wykorzystano metodę Runge-Kutty drugiego, trzeciego i czwartego rzędu:
Dokładność rozwiązania a dokładność punktów początkowych
Jak widać na powyższych wykresach, w przypadku metody Milne a (i ogólnie w przypadku wszystkich metod
wielokrokowych) obliczenie wartości początkowych ze zbyt małą dokładnością może uczynić całe rozwiązanie
bezużytecznym.
Metody "półanalityczne"
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 22
Rozważmy zagadnienie Cauchy-ego dla równania różniczkowego z warunkiem początkowym:
y f x,y ; y x0 y0
Metoda kolejnych przybliżeń polega na znalezieniu y(x) jako granicy ciągu funkcji yn(x) znajdowanych za pomocą
wzoru rekurencyjnego:
x
yn x y0 f x ,yn 1 x dx
x0
Udowodniono, że jeśli w pewnym obszarze R (określonym jako: x  x0 a oraz y  y0 b) spełniony jest warunek:
f x,y1 f x,y2 N y1 y2
gdzie N jest jakąś stałą, to niezależnie od wyboru funkcji początkowej y0(x) kolejne przybliżenia yn(x) w pewnym
przedziale (x,x+h) dążą do poszukiwanego rozwiązania zagadnienia.
Błąd n-tego przybliżenia w tym przedziale można oszacować następująco:
n 1
x x0
n
µn y x yn x MN
n 1 !
gdzie M jest określone jako:
M max f x,y
x,y R
a h jest określone jako:
b
h min a,
M
Metodę tę można również łatwo zastosować do układów równań różniczkowych.
Przykład: Znalez
ć kilka kolejnych przybliżeń rozwiązania równania: y = x2 + y2 dla warunku początkowego y(0) = 0.
Zapiszmy równanie w postaci całkowej:
x
2 2
y x x y dx
0
Jako przybliżenie wyjściowe przyjmiemy y0(x) 0. Kilka kolejnych przybliżeń będzie mieć następującą postać:
y0 0
3
x
y1
3
3 7
x x
y2
3 63
3 7 11 15
x x 2x x
y3
3 63 2079 59535
3 7 11 15 19
x x 2x 13x 82x
y4
3 63 2079 218295 37328445
23 27 31
662x 4x x
10438212015 3341878155 109876902975
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 23
3 7 11 15 19
x x 2x 13x 46x
y5
3 63 2079 218295 12442815
23 27 31
7382x 428x 17853193x
39665205657 50998121559 52289648987481675
35 39
738067298x 10307579354x
60394544580541334625 26851414520508677374275
43 47
6813116x 89797289962x
667791107219128471425 407191917287292044684671125
51 55
1159084x 721012x
308193360994502430591375 15328564533673936679413125
59 63
8x x
21664518085681213656375 760594829864786522589375
Określimy błąd przybliżenia numer 3. Przyjmijmy a = 1 i b = 0,5 (mogą to być dowolne liczby jako że funkcja x2 + y2
jest określona i ciągła na całej płaszczyz
nie xy). Otrzymamy wówczas:
2 2
M max f x,y max x y 1,25
N max f x,y max 2y 1
skąd wynika, że h = 0,4. Dla przedziału (0;0,4) otrzymamy:
4
x
5
4
y x y3 x 1,25 13 x
96
4!
i dla całego przedziału:
5 0,43
m a x y x y3 x 0,00133
96
(0;0,4)
Obliczone w powyższy sposób błędy poszczególnych przybliżeń są następujące: n=2  0.1, n=3  0.01(3), n=4
 0.001(3), n=5  0.00010(6). Jak widać dla każdego kolejnego przybliżenia błąd maleje o rząd wielkości.
Zastosowany wzór na błąd często daje wynik zbyt ostrożny (tzn. obliczony błąd jest zbyt duży). W praktyce obliczenia
prowadzi się do chwili, gdy kolejne yn i yn 1 zgadzają się z wymaganą dokładnością.
Na rysunku przedstawiającym różnice kolejnych przybliżeń znacznie powiększono skalę osi y. Oś x obejmuje nieco
zbyt duży zakres  do obliczeń przyjęto zakres x-ów od  1 do +1. Jeśli chcemy aby rozwiązanie było dokładne dla
większego zakresu x musimy posłużyć się przybliżeniem wyższego rzędu. W sytuacji dostępności programów całkujących
metoda ta może być bardzo użyteczna.
Przykład: Zamiast przyjąć jako funkcję wyjściową funkcję tożsamościowo równą zeru (jak w poprzednim przykładzie)
można znalez
ć funkcję bardziej zbliżoną do poszukiwanego rozwiązania. Na przykład dla równania y = x + 0,1y2 z
warunkiem początkowym y(0) = 1 funkcją wyjściową może być:
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 24
y 0 y 0
2
y0 x y0 x x
1! 2!
Powyższy wzór uzyskaliśmy rozwijając y(x) wokół punktu x0=0 (wówczas x jest odległością od x0) w szereg Taylora
ograniczony do drugich pochodnych. Z warunków początkowych i z równania otrzymamy wówczas:
y0 1
2
y 0 0,1y0 0,1
y 0 1 0,2y0y0 1,02
zatem ostatecznie będzie to funkcja:
2
y0 x 1 0,1x 0,51x
Zagadnienia brzegowe
W przypadku zagadnień początkowych dysponowaliśmy pełną informacją o szukanej funkcji w jednym wybranym
punkcie początkowym (tzn. znaliśmy wartość funkcji i wartości odpowiedniej liczby jej pochodnych  zależnie od
stopnia rozwiązywanego równania). W przypadku zagadnień brzegowych taką informacją nie dysponujemy, jest ona
 rozłożona pomiędzy dwa punkty ograniczające interesujący nas obszar, w którym chcemy znalez
ć rozwiązanie.
Wprawdzie rozwiązywanie zagadnień początkowych jest w pewnym sensie łatwiejsze (możliwe jest zastosowanie dużo
prostszych metod) to jednak kolejne obliczane wartości funkcji są obarczone coraz to większymi błędami, w krańcowych
przypadkach może to powodować konieczność porzucenia prostych metod rozwiązywania. W przypadku zagadnień
brzegowych ograniczenia narzucone na oba krańce przedziału powodują bardziej równomierne rozłożenie błędów
znalezionych rozwiązań.
Metoda strzałów
Metoda ta jest logiczną konsekwencją chęci zastosowania do zagadnień brzegowych metod opracowanych dla
zagadnień początkowych (a więc w zasadzie nie nadających się do takiego zastosowania). Czasami nazywana jest metodą
prób i błędów.
Załóżmy, że interesuje nas zagadnienie brzegowe dla równania drugiego rzędu. W przypadku zagadnienia
początkowego oczekiwalibyśmy wartości funkcji i jej pochodnej w punkcie początkowym, natomiast w przypadku
zagadnienia brzegowego będziemy dysponować wartościami funkcji na obu krańcach przedziału (oznaczmy je
odpowiednio przez a i b oraz przez f(a) i f(b)) natomiast nie będzie dana żadna wartość pochodnej szukanej funkcji..
Jeśli odgadniemy (obojętne czy biorąc wartość z sufitu czy
też szacując na podstawie jakichś sensownych przesłanek)
wartość pochodnej na jednym z krańców (oznaczmy ją przez
c1), to będziemy mogli zastosować którąś z wcześniej
opisanych metod i za jej pomocą obliczyć wszystkie wartości
szukanej funkcji aż do drugiego brzegu włącznie (wartość tę
oznaczymy przez f1). Nie oczekujemy oczywiście, że obliczona
w ten sposób wartość funkcji w punkcie b będzie zgadzać się
z wynikającą z warunków zadania wartością brzegową. Jednak
uzyskamy w ten sposób parę wartości: założona wartość
pochodnej i otrzymana wartość funkcji na drugim brzegu.
Teraz możemy przyjąć inną wartość pochodnej na brzegu (c2)
i znów przeprowadzić obliczenia (otrzymując f2).
Jeśli będziemy mieć szczęście jedna z otrzymanych
wartości będzie za duża a druga za mała. Wtedy wartość c3 wyznaczymy choćby za pomocą interpolacji liniowej:
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 25
c2 c1 f1 f b
c3 c1
f1 f2
(przy założeniu, że f1 było za duże a f2 za małe).
Pierwszy przykład nie jest właściwie zagadnieniem brzegowym  rozwiązywane jest równanie pierwszego rzędu,
jednak podano warunek (wartość funkcji szukanej) nie dla początku lecz dla końca przedziału. Można byłoby oczywiście
rozwiązać to zagadnienie którąkolwiek z poznanych już metod stosując krok o ujemnej wartości jednak i w takim
przypadku można zastosować metodę strzałów.
Kolejne przykłady dotyczą równania drugiego rzędu z podanymi wartościami funkcji szukanej na obu krańcach
przedziału (jest to już więc zagadnienie brzegowe). Jak widać całkiem niezłą dokładność rozwiązania można uzyskać
nawet bez nadmiernej dbałości o dokładność poszukiwanej wartości początkowej pierwszej pochodnej.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 26
Metody różnicowe
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 27
Metody różnicowe w odniesieniu do zagadnień brzegowych dla równań różniczkowych zwyczajnych są takie same
jak w przypadku równań różniczkowych cząstkowych, rozwiązania są jednak prostsze, bo mamy do czynienia tylko z
jedną zmienną niezależną. Zastosowane mogą być praktycznie wszystkie wzory podane dla równań różniczkowych
czÄ…stkowych.
Metoda passage
Dla równania liniowego drugiego rzędu możliwe jest wyprowadzenie wzorów z uwzględnieniem dowolnych warunków
brzegowych (trzeciego rodzaju  zawierających i wartość funkcji i wartość pochodnej na brzegu przedziału). Załóżmy,
że rozwiązywany problem ma postać:
y p(x)y q(x)y f(x)
Ä…0y(a) Ä…1y (a) A
²0y(b) ²1y (b) B
gdzie p(x), q(x) i f(x) sÄ… znanymi funkcjami, a i b sÄ… wartoÅ›ciami granicznymi zmiennej x zaÅ› Ä…0, Ä…1, ²0, ²1, A i B sÄ…
znanymi stałymi. Zakładamy dalej, że odcinek a b został podzielony na N kroków o długości "x, przy czym węzeł nr 0
odpowiada lewej granicy, a węzeł nr N  granicy prawej.
Wówczas oznaczając:
2
2qi "x 4 2 pi "x
mi ; ki
2 pi "x 2 pi "x
możemy przeprowadzić obliczenia jak następuje. Najpierw obliczamy stałe c1 i d1:
Ä…1 Ä…0"x
c1
m1 Ä…1 Ä…0"x k1Ä…1
2
2f1"x
d1 k1 A"x
2 p1"x Ä…1 Ä…0"x
następnie wyznaczamy pozostałe stałe aż do cN i dN:
1
ci
mi kici 1
2
2fi "x
di kici 1di 1
2 pi "x
Na podstawie cN i dN obliczamy yN:
2B"x ²1 dN cN 1dN 1
yN
1
2²0"x ²1 cN 1
cN
po czym kolejne wartości yi (to, że wzory przewidują możliwość obliczenia wartości yo wynika z faktu, że warunki
brzegowe mogą nie zawierać wartości funkcji na brzegu  gdy nie są to warunki pierwszego rodzaju):
yi ci di yi 1
Ä…1y1 A"x
y0
Ä…1 Ä…0"x
Wyprowadzenie: Rozwiązywane równanie ma postać:
y p x y q x y f x
a warunki brzegowe:
Ä…0y a Ä…1y a A
²0y b ²1y b B
W postaci różnicowej:
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 28
yi 1 2yi yi 1 yi 1 yi 1
pi qiyi fi
2
2h
h
y1 y0
Ä…0y0 Ä…1 A
h
yN yN
²0yN ²1 1 1 B
2h
Warto zauważyć, że do zapisania prawego warunku brzegowego (dla x=b) wykorzystano różnicę centralną  zapis
różnicowy nie jest więc symetryczny względem x i należy oczekiwać, że wyniki mogą być różne w zależności od tego od
którego końca przedziału zaczniemy rozwiązywać.
Po przekształceniu otrzymamy:
yi 1 miyi kiyi 1 i
gdzie:
2 2
2qih 4 2 hpi 2h fi
mi ; ki ;
i
2 hpi 2 hpi 2 hpi
Jest to oczywiście układ równań opisany macierzą trójdiagonalną. Załóżmy, że dla i = 1,2,...,n 1 możliwe jest zapisanie
rozwiązań w postaci:
yi ci di yi 1
Jeśli jest to prawdą, to możemy podstawić taki wzór na yi 1 do różnicowej postaci rozwiązywanego równania:
yi 1 ci 1 di 1 yi
yi 1 mi yi ki ci 1 di 1 yi i
co po przekształceniach daje:
1
yi ki ci di yi
mi ki ci 1 i 1 1 1
czyli:
1
ci ; di ki ci 1di 1
i
mi ki ci 1
Dla i = 1 otrzymamy:
y2 m1y1 k1y0 1
Z warunku brzegowego dla x = a mamy:
y1 y0 Ah Ä…1y1
Ä…0y0 Ä…1 A y0
h Ä…0h Ä…1
po podstawieniu:
Ah Ä…1y1
y2 m1y1 k1
Ä…0h Ä…1 1
i po przekształceniu:
Ä…1k1 k1Ah
m1 y1 y2
1
Ä…oh Ä…1 Ä…0h Ä…1
Ä…0h Ä…1 k1Ah
y1 y2
1
m1 Ä…0h Ä…1 Ä…1k1 Ä…0h Ä…1
co daje nam wartości c1 i d1. Na podstawie wzorów na ci i di możemy obliczyć wszystkie kolejne wartości współczynników
c i d. Następnie rozważymy układ równań (który umożliwi nam wyeliminowanie z obliczeń wartości yN+1):
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 29
yN cN dN yN 1
yN 1 cN 1 dN 1 yN
yN yN
²0 yN ²1 1 1 B
2h
po przekształceniach otrzymamy:
yN
yN 1 dN
cN
²1 yN
²0 yN dN cN 1dN 1 cN 1yN B
2h cN
i dalej:
1
2h²0 ²1 cN 1 yN 2hB ²1 cN 1dN 1 dN
cN
Ostatecznie wartość yN wyniesie:
2hB ²1 cN 1dN 1 dN
yN
1
2h²0 ²1 cN 1
cN
Po obliczeniu yN możemy  cofając się obliczyć wszystkie pozostałe wartości y, aż do y1.
Lewy warunek brzegowy zapisany w postaci różnicowej możemy przekształcić do postaci umożliwiającej obliczenie
y0:
Ah Ä…1y1
y0
Ä…1 Ä…0h
Metoda różnicowa dla równań nieliniowych drugiego rzędu
Rozważmy teraz zagadnienie brzegowe dla nieliniowego równania drugiego rzędu:
y f x,y,y
z warunkami brzegowymi:
Ä…0 y a Ä…1y a A
²0 y b ²1y b B
Równanie różniczkowe zastąpimy następującym równaniem różnicowym:
yi 1 2yi yi 1 yi 1 yi 1
f xi,yi,
2
2h
h
y1 y0 yN yN
Ä…0 y0 Ä…1 A ; ²0 yN ²1 1 B
h h
Jest to równanie nieliniowe, rozwiązujemy je metodą iteracyjną. Oznaczmy:
y1 y0
“0 y Ä…0y0 Ä…1
h
yN yN
“N y ²0yN ²1 1
h
iteracje będą mieć postać (indeks górny (r) oznacza numer iteracji):
yi(r 1) 2yi(r 1) yi(r 1) yi(r ) yi(r )
1 1 1 1
f xi,yi(r ), ; i 1, ,N 1
2
2h
h
“0 yi(r 1) A (i 0) ; “N yi(r 1) B (i N )
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 RozwiÄ…zywanie ODE 30
W ramach każdej z iteracji musimy rozwiązać układ równań algebraicznych. Rozwiązanie to można przedstawić w
postaci:
N 1
h i (r)
2
yi(r 1) A²0 b a A ²1 Ä…1B Ä…0B A²0 h gji f
j
" "
j 1
gdzie gji i " obliczane sÄ… jako:
1
" Ä…0²0 b a Ä…0²1 Ä…1²0
h
Ä…1 ²1
1
j Ä…0 i ²0 ²0N j i
" h h
gji
Ä…1 ²1
1
i Ä…0 j ²0 ²0N j>i
" h h
Warto zauważyć, że we wzorze iteracyjnym tylko wartość fj zależy od numeru iteracji.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej