Fizyka laboratorium 3
1.1. Ró żnice skończone
Przyjrzyjmy sie równaniu różniczkowemu zwyczajnemu pierwszego rzedu, postaci:
,
,
dy(t)
f (y, t) =
(1.1)
dt
Równanie (1.1) opisuje zmiane funkcji y przy zmianie parametru t, Jeżeli parametr
,
t bedzie oznacza l czas, równanie to opisuje wyrażona przez funkcje f zmiane funkcji
,
,
,
,
y w czasie. Rozwiazaniem analitycznym jest tego równania jest funkcja y(t). która
,
podstawiona do równania (1.1) przekszta lca je w tożsamość.
Zamieniajac odpowiednio przyrost dy na ∆y oraz dt na ∆t otrzymujemy:
,
∆y
dy
lim
=
(1.2)
∆y→0,∆t→0
∆t
dt
Równanie (1.1) w rachunku różnic skończonych zapisać możemy w postaci:
∆y
f (y, t) =
(1.3)
∆t
Mnożac obie strony przez ∆t otrzymujemy wyrażenie na zmiane funkcji y w prze-
,
,
dziale czasu ∆t:
∆y = f (y, t)∆t
(1.4)
Znajac warunki poczatkowe (poczatkowa wartość y) możemy obliczyć nowa wartość
,
,
,
,
y po up lywie czasu ∆t jako:
∆y
= f (y, t)∆t
(1.5)
yn+1 = yn + f (yn, yn)∆t
1
1.2. Ca lkowanie równań różniczkowych 1.2. Ca lkowanie równań ró żniczkowych Równanie (1.1) przedstawia funkcje y zależna od czasu t. Równanie to można rozwia-
,
,
,
zać poprzez ca lkowanie po parametrze t. Ponieważ bedziemy rozwiazywali równanie
,
,
dla kolejnych kroków czasowych wprowadzimy nastepujace oznaczenia: y
,
,
n+1 – war-
tość y w n + 1 kroku czasowym, yn – wartość y w n-tym kroku czasowym: Z
tn+1
yn+1 = yn +
f (y, t)dt
(1.6)
tn
Równanie to opisuje przyrost funkcji y w przedziale czasu od tn do tn+1. Ca la idea metody numerycznej skupia sie na przybliżeniu różnicowym ca lki po prawej stro-
,
nie równania. Dok ladność, z jaka wyznaczymy wartość ca lki, ma znaczenie dla sta-
,
bilności rozwiazania numerycznego równania (1.1). Należy pamietać, że otrzymane
,
,
rozwiazanie jest tylko przybliżeniem rozwiazania dok ladnego.
,
,
1.3. Schemat ró żnicowy Eulera
Najprostszym sposobem na przybliżenie ca lki z równania (1.6) jest ca lkowanie metoda prostokatów.
,
,
f(y,t)
fn
t
tn
tn+1
Rysunek 1.1. metoda prostokatów
,
Stosujac te metode otrzymamy najprostszy schemat różnicowy, zwany metoda Eu-
,
,
,
,
lera:
yn+1 = yn + f (y, t)h + O(h2)
(1.7)
gdzie h ↔ ∆t
Możemy wyprowadzić schemat metody Eulera siegajac do definicji pochodnej:
,
,
0
y(t + h) − y(t)
y (t) = lim
(1.8)
h→0
h
możemy wiec rezygnujac z granicy możemy zapisać przybliżenie:
,
,
0
y(t + h) − y(t)
y (t) ≈
(1.9)
h
2
1.4. Metoda Midpoint czyli po przekszta lceniu
0
y(t + h) − y(t) ≈ hy (t)
(1.10)
czyli:
0
y(t + h) ≈ y(t) + hy (t)
(1.11)
0
Przechodzac do wprowadzonej notacji y (t) = f (y, t) możemy zapisać schemat w
,
postaci:
y(t + h) ≈ y(t) = hf (y, t)
(1.12)
yn+1
styczna
f ’(x n )
f(x n+1 )
f(x n )
h=x n+1- x n
x n
x n+1
Rysunek 1.2. metoda Eulera
1.3.1. Rzad dok ladności metody
,
Dla sprawdzenia rzedu dok ladności metody dokonamy rozwiniecia funkcji y w szereg
,
,
Taylora:
∞
00
X y(n)
0
y (t)
y(n)
y(t + h) =
hn = y(t) + y (t)h +
h2 + . . . +
hn
(1.13)
n!
2!
n!
n=0
Widać, że metoda Eulera obcina rozwiniecie po pierwszej pochodnej czyli: y(t+h) ≡
,
0
yn+1, y(t) ≡ yn, y (t)h ≡ f (y, t)h Obciecie rozwiazania do wyrazu h w potedze 1 w lacznie oznacza, że metoda różni-
,
,
,
,
cowa ma dok ladność pierwszego rzedu wzgledem rozwiniecia w szereg Taylora. B lad
,
,
,
,
obciecia wynosi O(h2)
,
1.4. Metoda Midpoint
Nazwa metody wynika z faktu, że tn + h/2 jest punktem środkowym pomiedzy
,
tn w którym wartość funkcji y(t) jest znana, oraz punktem tn+1 w którym szukamy 3
1.4. Metoda Midpoint wartości funkcji y(t). Jest to metoda zaliczana do metod Runge-Kutta (inaczej zwana metoda Runge-Kutta drugiego rzedu).
,
Metoda ta bazuje na metodzie Eulera wykorzystujac kolejna pochodna przy rozwi-
,
,
,
nieciu w szereg Taylora.
,
0
h
y(t + h) − y(t)
y (t +
) ≈
(1.14)
2
h
nastepnie:
,
0
h
y(t + h) ≈ y(t) + hy (t +
)
(1.15)
2
co możemy zapisać jako:
h
h
y(t + h) ≈ y(t) + hf
y(t +
), t +
(1.16)
2
2
Nie możemy skorzystać z równania (1.16) ponieważ nie znamy wartości y w punkcie h
t +
. Rozwiazaniem jest rozwiniecie w szereg Taylora i przybliżenie tej wartości.
2
,
,
h
h 0
h
y
t +
≈ y(t) +
y (t) = y(t) +
f (y(t), t)
(1.17)
2
2
2
Co po podstawieniu do (1.16) daje:
h
h
y(t + h) ≈ y(t) + hf
y(t) +
f (y(t), t) , t +
(1.18)
2
2
Ogólnie schemat metody MidPoint możemy zapisać jako: k1 = f (yn, tn)
(1.19)
h
h
k2 = f
yn +
k1, tn +
(1.20)
2
2
yn+1 = yn + hk2 + O(h3)
(1.21)
4
1.5. Metoda Runge–Kutta Przybli•ona styczna
wynik Euler yn+1
w punkcie t +h/2
n
k = f ( y , t ) 1
n
n
h
y +
k
n
1
2
wynik MidPoint
yn 1
+
yn
h/2
h/2
t
t
n
n 1
+
Rysunek 1.3.
metoda MidPoint
1.5. Metoda Runge–Kutta
Ogólnie metode możemy zapisać w postaci:
,
X
yn+1 = yn + h
amkm
(1.22)
m
gdzie: am – wspó lczynniki wagowe, km styczne obliczone w kolejnych krokach tn ≤
x ≤ tn+1. Wagi spe lniaja warunek:
,
X am = 1
(1.23)
m
1.5.1. Metoda RK4
Schematycznie metode możemy zapisać w nastepujacy sposób:
,
,
,
k1 = f (yn, tn)
h
h
k2 = f
yn +
k1, tn +
2
2
h
h
k3 = f
yn +
k2, tn +
2
2
k4 = f (yn + hk3, tn + h)
co zgodnie ze wzorem (1.22) daje: h
yn+1 = yn +
(k1 + 2k2 + 2k3 + k4)
(1.24)
6
5
1.5. Metoda Runge–Kutta wynik Euler yn+1
Przybli•ona styczna
w punkcie t +h/2
n
k = f ( y , t ) 1
n
n
h
h
h
y +
k
n
1
=
+
+
2
k
f ( y
k , x
2
n
2 1
n
2
2
h
h
y
k
n 1
+
= f ( y + k , x +
3
n
2 2
n
2
3
yn
=
+
+
1
k
f ( y
hk , t
h)
4
n
n
3
h/2
h/2
t
t
n
n 1
+
Rysunek 1.4. metoda RK4
6