Rozdzia l 1

Fizyka laboratorium 3

1.1. R´

o ˙znice sko´

nczone

Przyjrzyjmy si

,

e r´

ownaniu r´

o˙zniczkowemu zwyczajnemu pierwszego rz

,

edu, postaci:

f (y, t) =

dy(t)

dt

(1.1)

R´

ownanie (

1.1

) opisuje zmian

,

e funkcji y przy zmianie parametru t, Je˙zeli parametr

t b

,

edzie oznacza l czas, r´

ownanie to opisuje wyra˙zon

,

a przez funkcj

,

e f zmian

,

e funkcji

y w czasie. Rozwi

,

azaniem analitycznym jest tego r´

ownania jest funkcja y(t). kt´

ora

podstawiona do r´

ownania (

1.1

) przekszta lca je w to˙zsamo´s´

c.

Zamieniaj

,

ac odpowiednio przyrost dy na ∆y oraz dt na ∆t otrzymujemy:

lim

∆y→0,∆t→0

 ∆y

∆t



=

dy

dt

(1.2)

R´

ownanie (

1.1

) w rachunku r´

o˙znic sko´

nczonych zapisa´

c mo˙zemy w postaci:

f (y, t) =

∆y

∆t

(1.3)

Mno˙z

,

ac obie strony przez ∆t otrzymujemy wyra˙zenie na zmian

,

e funkcji y w prze-

dziale czasu ∆t:

∆y = f (y, t)∆t

(1.4)

Znaj

,

ac warunki pocz

,

atkowe (pocz

,

atkowa warto´s´

c y) mo˙zemy obliczy´

c now

,

a warto´s´

c

y po up lywie czasu ∆t jako:

 ∆y

= f (y, t)∆t

y

n+1

= y

n

+ f (y

n

, y

n

)∆t

(1.5)

1

1.2. Ca lkowanie r´

owna´

n r´

o˙zniczkowych

1.2. Ca lkowanie r´

owna´

n r´

o ˙zniczkowych

R´

ownanie (

1.1

) przedstawia funkcj

,

e y zale˙zn

,

a od czasu t. R´

ownanie to mo˙zna rozwi

,

a-

za´

c poprzez ca lkowanie po parametrze t. Poniewa˙z b

,

edziemy rozwi

,

azywali r´

ownanie

dla kolejnych krok´

ow czasowych wprowadzimy nast

,

epuj

,

ace oznaczenia: y

n+1

– war-

to´s´

c y w n + 1 kroku czasowym, y

n

– warto´s´

c y w n-tym kroku czasowym:

y

n+1

= y

n

+

Z

t

n+1

t

n

f (y, t)dt

(1.6)

R´

ownanie to opisuje przyrost funkcji y w przedziale czasu od t

n

do t

n+1

. Ca la idea

metody numerycznej skupia si

,

e na przybli˙zeniu r´

o˙znicowym ca lki po prawej stro-

nie r´

ownania. Dok ladno´s´

c, z jak

,

a wyznaczymy warto´s´

c ca lki, ma znaczenie dla sta-

bilno´sci rozwi

,

azania numerycznego r´

ownania (

1.1

). Nale˙zy pami

,

eta´

c, ˙ze otrzymane

rozwi

,

azanie jest tylko przybli˙zeniem rozwi

,

azania dok ladnego.

1.3. Schemat r´

o ˙znicowy Eulera

Najprostszym sposobem na przybli˙zenie ca lki z r´

ownania (

1.6

) jest ca lkowanie me-

tod

,

a prostok

,

at´

ow.

f(y,t)

t

n

t

n+1

f

n

t

Rysunek 1.1. metoda prostok

,

at´

ow

Stosuj

,

ac t

,

e metod

,

e otrzymamy najprostszy schemat r´

o˙znicowy, zwany metod

,

a Eu-

lera:

y

n+1

= y

n

+ f (y, t)h + O(h

2

)

(1.7)

gdzie h ↔ ∆t

Mo˙zemy wyprowadzi´

c schemat metody Eulera si

,

egaj

,

ac do definicji pochodnej:

y

0

(t) = lim

h→0

y(t + h) − y(t)

h

(1.8)

mo˙zemy wi

,

ec rezygnuj

,

ac z granicy mo˙zemy zapisa´

c przybli˙zenie:

y

0

(t) ≈

y(t + h) − y(t)

h

(1.9)

2

1.4. Metoda Midpoint

czyli po przekszta lceniu

y(t + h) − y(t) ≈ hy

0

(t)

(1.10)

czyli:

y(t + h) ≈ y(t) + hy

0

(t)

(1.11)

Przechodz

,

ac do wprowadzonej notacji y

0

(t) = f (y, t) mo˙zemy zapisa´

c schemat w

postaci:

y(t + h) ≈ y(t) = hf (y, t)

(1.12)

x

n

x

n+1

f(x

n

)

f(x

n+1

)

y

n+1

styczna

f

’

(x

n

)

h=x

n+1

-x

n

Rysunek 1.2. metoda Eulera

1.3.1. Rz

,

ad dok ladno´

sci metody

Dla sprawdzenia rz

,

edu dok ladno´sci metody dokonamy rozwini

,

ecia funkcji y w szereg

Taylora:

y(t + h) =

∞

X

n=0

y

(n)

n!

h

n

= y(t) + y

0

(t)h +

y

00

(t)

2!

h

2

+ . . . +

y

(n)

n!

h

n

(1.13)

Wida´

c, ˙ze metoda Eulera obcina rozwini

,

ecie po pierwszej pochodnej czyli: y(t+h) ≡

y

n+1

, y(t) ≡ y

n

, y

0

(t)h ≡ f (y, t)h

Obci

,

ecie rozwi

,

azania do wyrazu h w pot

,

edze 1 w l

,

acznie oznacza, ˙ze metoda r´

o˙zni-

cowa ma dok ladno´s´

c pierwszego rz

,

edu wzgl

,

edem rozwini

,

ecia w szereg Taylora. B l

,

ad

obci

,

ecia wynosi O(h

2

)

1.4. Metoda Midpoint

Nazwa metody wynika z faktu, ˙ze t

n

+ h/2 jest punktem ´srodkowym pomi

,

edzy

t

n

w kt´

orym warto´s´

c funkcji y(t) jest znana, oraz punktem t

n+1

w kt´

orym szukamy

3

1.4. Metoda Midpoint

warto´sci funkcji y(t). Jest to metoda zaliczana do metod Runge-Kutta (inaczej zwana
metoda Runge-Kutta drugiego rz

,

edu).

Metoda ta bazuje na metodzie Eulera wykorzystuj

,

ac kolejn

,

a pochodn

,

a przy rozwi-

ni

,

eciu w szereg Taylora.

y

0

(t +

h

2

) ≈

y(t + h) − y(t)

h

(1.14)

nast

,

epnie:

y(t + h) ≈ y(t) + hy

0

(t +

h

2

)

(1.15)

co mo˙zemy zapisa´

c jako:

y(t + h) ≈ y(t) + hf



y(t +

h

2

), t +

h

2



(1.16)

Nie mo˙zemy skorzysta´

c z r´

ownania (

1.16

) poniewa˙z nie znamy warto´sci y w punkcie

t +

h

2

. Rozwi

,

azaniem jest rozwini

,

ecie w szereg Taylora i przybli˙zenie tej warto´sci.

y



t +

h

2



≈ y(t) +

h

2

y

0

(t) = y(t) +

h

2

f (y(t), t)

(1.17)

Co po podstawieniu do (

1.16

) daje:

y(t + h) ≈ y(t) + hf



y(t) +

h

2

f (y(t), t) , t +

h

2



(1.18)

Og´

olnie schemat metody MidPoint mo˙zemy zapisa´

c jako:

k

1

= f (y

n

, t

n

)

(1.19)

k

2

= f



y

n

+

h

2

k

1

, t

n

+

h

2



(1.20)

y

n+1

= y

n

+ hk

2

+ O(h

3

)

(1.21)

4

1.5. Metoda Runge–Kutta

1

+

n

y

1

2

k

h

y

n

+

wynik Euler y

n+1

h/2

h/2

Przybli•ona styczna
w punkcie t

n

+h/2

wynik MidPoint

)

,

(

1

n

t

n

y

f

k

=

n

y

n

t

1

+

n

t

Rysunek 1.3. metoda MidPoint

1.5. Metoda Runge–Kutta

Og´

olnie metod

,

e mo˙zemy zapisa´

c w postaci:

y

n+1

= y

n

+ h

X

m

a

m

k

m

(1.22)

gdzie: a

m

– wsp´

o lczynniki wagowe, k

m

styczne obliczone w kolejnych krokach t

n

≤

x ≤ t

n+1

. Wagi spe lniaj

,

a warunek:

X

m

a

m

= 1

(1.23)

1.5.1. Metoda RK4

Schematycznie metod

,

e mo˙zemy zapisa´

c w nast

,

epuj

,

acy spos´

ob:

k

1

= f (y

n

, t

n

)

k

2

= f



y

n

+

h

2

k

1

, t

n

+

h

2



k

3

= f



y

n

+

h

2

k

2

, t

n

+

h

2



k

4

= f (y

n

+ hk

3

, t

n

+ h)

co zgodnie ze wzorem (

1.22

) daje:

y

n+1

= y

n

+

h

6

(k

1

+ 2k

2

+ 2k

3

+ k

4

)

(1.24)

5

1.5. Metoda Runge–Kutta

h/2

h/2

1

2

3

1

2

k

h

y

n

+

)

,

3

(

4

2

,

2

2

(

3

2

,

1

2

(

2

)

,

(

1

h

n

t

hk

n

y

f

k

h

n

x

k

h

n

y

f

k

h

n

x

k

h

n

y

f

k

n

t

n

y

f

k

+

+

=

+

+

=

+

+

=

=

n

t

1

+

n

t

n

y

1

+

n

y

Przybli•ona styczna
w punkcie t

n

+h/2

wynik Euler y

n+1

Rysunek 1.4. metoda RK4

6