Równania różniczkowe zwyczajne

(1)

Równania różniczkowe zwyczajne

W najprostszym przypadku poszukujemy funkcji jednej zmiennej rzeczywistej x w postaci:

( )

x y

y= , (1)

której pochodna y'

( )

x ma spełniać równanie dane jako:

( )

x f

(

x y

( )

x

)

y' = , , (2)

lub w skrócie:

( )

x y f

y'= , . (3)

W ogólności tak postawiony problem ma nieskończenie wiele rozwiązań, spomiędzy których wy- bieramy rozwiązanie szczególne, spełniające dodatkowy warunek, zwany warunkiem początkowym – poszukujemy takiej funkcji (1), która dla danego x spełnia: ₀

( )

x₀ y₀

y = ; (4)

Zadanie (2), (4) można uogólnić na przypadki:

• układu n równań różniczkowych zwyczajnych:

( ) ( ( ) ( ) ( ) ( ) )

⎪⎪

⎩

⎪⎪

⎨

⎧

=

x y x y x y x y x f x y

n i

n n

n i

,..., ,...,

, ,

,..., ,...,

, ,

,..., ,...,

, ,

2 1 '

2 1 2 '

2

2 1 1 '

1

M , (5)

n nieznanych funkcji rzeczywistych y_i

( )

x, i=1,2,...,i,...,n zmiennej rzeczywistej x z warunkiem początkowym:

( ) ( )

( )

₀ ₀

20 0 2

10 0 1

n

n x y

y

y x y

=

M ; (6)

• równania różniczkowego rzędu wyższego niż pierwszy:

( )

^x ^f

(

^x ^y

( ) ( ) ( )

^x ^y ^x ^y ^x ^y⁽ ⁾

( )

^x

)

y^m = , , ^' , ^'' ,..., ^m⁻¹ , (7)

jednej funkcji rzeczywistej (1) z warunkiem początkowym:

( )

x₀ = y₀ ,y^'

( )

x₀ = y₀^' ,y^''

( )

x₀ =y₀^'',...,y⁽^m⁻¹⁾

( )

x₀ = y₀⁽^m⁻¹⁾

y . (8)

Przypadek drugi z rozważanych powyżej zawsze można sprowadzić do przypadku pierwszego po- przez ciąg podstawień:

( ) ( ) ( ) ( ) ( )

⁽ ⁾

( )

( ) ( ) ( ) ( ) ( ) ( )

( )

x f

(

x z

( ) ( )

x z x z

( )

x

)

z

x z x z

x y x z

m m

m

m ^' , ₁ , ₂ ,...,

' 1

3 '

2

2 '

1

1 ' 2

1

=

⇒

=

− −

M M . (9)

(2)

Politechnika Krakowska

Ze względu na sposób postępowania numerycznego, metody rozwiązania problemu (2) z warunkiem (4) można podzielić na dwie kategorie:

• metody jednokrokowe (znane również pod nazwą samostartujących), należą tu na przykład metoda Eulera, bądź metody typu Rungego – Kutty,

• metody wielokrokowe (znane także jako metody typu predyktor – korektor), do których na- leżą np. metody Adamsa – Bashforda (predyktor) i Adamsa – Moultona (korektor).

Metody jednokrokowe można zapisać przy pomocy następujących wzorów:

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛ + ⋅α + ⋅ β ⋅

=

∑

= s

j ij j

n i n

i f x h y h k

k

1

, , (10)

0 ,

1

1= + ⋅

∑

⋅ ≥

+ = i

s

i i i

n

n y h w k w

y , (11)

wynikających z rozwinięcia poszukiwanej funkcji (1) w otoczeniu punktu y w szereg Taylora. _n Konkretne wartości współczynników w , _i α i _i β w przypadku metod jawnych, dla których musi _ij być spełniony warunek s≤ we wzorze (10), przedstawiają się następująco: i

• dla metody Eulera w₁=1, α₁=0, β₁₁=0, czyli:

( )

1 1

1 ,

k h y y

y x f k

n n

⋅ +

=

+

; (12)

• dla metody Rungego – Kutty II rzędu w₁=0, w₂ =1, α₁ =0, α₂ = ₂¹, β₁₁=0 i β₂₂= ₂¹, czyli:

( )

(

₁ ₂

)

1

1 21 21

2

1

1 0 ,

,

k k h y y

k h y h x f k

y x f k

n n

⋅ +

⋅

⋅ +

=

⋅

⋅ +

=

+

; (13)

• dla metody Rungego – Kutty IV rzędu w₁ =₆¹, w₂ =¹₃, w₃ =₃¹, w₄= ₆¹, α₁=0, α₂ = ₂¹,

21 3 =

α , α₄ =1, β₁₁=0, β₂₂=¹₂, β₃₃= ₂¹ i β₄₄=1 czyli:

( )

(

¹₆ ₁ ₃¹ ₂ ₃¹ ₃ ₆¹ ₄

)

1

3 4

2 2 2 1

3 1

2 1 2 1

2 1

1

, , , ,

k k k k h y y

k h y h x f k

y x f k

n n

⋅ +

⋅

⋅ +

=

⋅ + +

=

⋅

⋅ +

=

⋅

⋅ +

=

+

. (14)

We wzorach (13) i (14) wszystkie współczynniki β , dla których _ij i≠ , są równe zeru. j

Metoda Eulera jest metodą rzędu pierwszego, czyli różnica pomiędzy ścisłą y_n^ść₊₁ a przybliżoną y _n₊₁ wartością rozwiązania w kolejnym punkcie x zmienia się z pierwszą potęgą _n₊₁ h . W metodzie Rungego – Kutty II rzędu różnica ta zmienia się z drugą, w metodzie IV rzędu z czwartą potęgą h . Dla ilustracji toku postępowania przy rozwiązywaniu problemu początkowego każdą z wymienio- nych powyżej metod, zastosujmy je do wykonania trzech kolejnych kroków w poszukiwaniu numerycznego rozwiązania zadania:

(3)

1 5 3

+

−

= ⋅ t

z dt

dz , (15)

przy warunku początkowym:

( )

2,0

, 0 ,

2 ₀ ₀

0 = z =z t =

t , (16)

z krokiem ∆t=h= ₂¹ . Problem ten ma rozwiązanie ścisłe w postaci funkcji z

( )

t danej wzorem:

( ) ( )

3 5 81

1³ + +

= t t

z^ść . (17)

Metoda Eulera Iteracja pierwsza:

( )

⁰^,³³³³³³

1 000 , 2

5 000000 ,

2 3 1

5 , 3

0 0 0

0

0 =

+

−

= ⋅ +

+

= ⋅

= t

z z t f

f , (18)

166667z₁=z₀+h⋅f₀=2,000000+0,500⋅0,333333=2, . (19) Iteracja druga:

( )

0,428571

1 500 , 2

5 166667 ,

2 3 1

5 , 3

1 1 1 1

1 =

+

−

= ⋅ +

+

= ⋅

= t

z z t f

f , (20)

380952 ,

2 428571 ,

0 500 , 0 166667 ,

1 2

1

2 =z +h⋅ f = + ⋅ =

z . (21)

Iteracja trzecia:

( )

0,535714

1 000 , 3

5 380952 ,

2 3 1

5 , 3

2 2 2

2

2 =

+

−

= ⋅ +

+

= ⋅

= t

z z t f

f , (22)

648810z₃=z₂+h⋅f₂ =2,380952+0,500⋅0,535714=2, . (23) Metoda Rungego – Kutty II rzędu

Iteracja pierwsza:

{ }

( )

{ }

(

{ }

) (

^{{ }}

)

( ) 0,384615

1 5 , 3

333333 ,

1 0 5 , 3

1 500 , 0 5 , 0 000 , 2

5 333333 , 0 500 , 0 5 , 0 000000 , 2 3 21

0

1 2 1 0 1 1

1 12 0 21 0 1 2

1 000 , 2

5 000000 , 2 3 0

0 0 0 1

1

= + =

⋅ +

−

⋅

⋅ +

= ⋅

⋅

⋅ +

=

= + =

−

= ⋅

=

+

⋅ +

−

⋅

⋅ +

−

⋅

h t

k h k z

h z h t f k

t z z t f k

, (24)

{ }¹ 2,000000 0,500 384615 2,192308

2 0

1=z +h⋅k = + ⋅ =

z . (25)

Iteracja druga:

{ }

( )

{ }

(

{ }

) (

^{{ }}

)

₍ ₎ ₀_,₅₁₀₆₂₃

1 5 , 3

450549 ,

1 0 5 , 3

1 500 , 0 5 , 0 500 , 2

5 450549 , 0 500 , 0 5 , 0 192308 , 2 3 21

1

2 1 12 2 1

1 21 1 21 1 2 2

1 500 , 2

5 192308 , 2 3 1

1 1 1 2 1

= + =

⋅ +

−

⋅

⋅ +

= ⋅

⋅

⋅ +

=

= + =

−

= ⋅

=

+

⋅ +

−

⋅

⋅ +

−

⋅

h t

k h k z

h z h t f k

t z z t f k

, (26)

{ }² 2,192308 0,500 0,510623 2,447619

2 1

2 =z +h⋅k = + ⋅ =

z . (27)

(4)

Iteracja trzecia:

{ }

( )

{ }

(

{ }

) (

^{{ }}

)

( ) 0,654622

1 5 , 3

585714 ,

1 0 5 , 3

1 500 , 0 5 , 0 000 , 3

5 585714 , 0 500 , 0 5 , 0 447619 , 2 3 12

2

3 1 21 3 2

2 1 2 1 21 2 3

2

1 000 , 3

5 447619 , 2 3 2

2 2 2 3

1

= + =

⋅ +

−

⋅

⋅ +

= ⋅

⋅

⋅ +

=

= + =

−

= ⋅

=

+

⋅ +

−

⋅

⋅ +

−

⋅

h t

k h k z

h z h t f k

t z z t f k

, (28)

{ }³ 2,447619 0,500 0,654622 2,774930

2 2

3 =z +h⋅k = + ⋅ =

z . (29)

Metoda Rungego – Kutty IV rzędu Iteracja pierwsza:

{ }

( )

{ }

(

{ }

) (

^{{ }}

)

( )

{ }

(

{ }

) (

^{{ }}

)

( )

{ }

(

{ }

) (

^{{ }}

)

( ) 0,455621

1 5 , 3

396450 ,

1 0 5 , 3

384615 ,

1 0 5 , 3

333333 ,

1 0 5 , 3

1 500 , 0 000 , 2

5 396450 , 0 500 , 0 000000 , 2 3 0

1 3 1 0

3 0 0 1

4

1 500 , 0 5 , 0 000 , 2

5 384615 , 0 500 , 0 5 , 0 000000 , 2 3 21

0

1 2 21 1 0

2 12 0 21 0 1 3

1 500 , 0 5 , 0 000 , 2

5 333333 , 0 500 , 0 5 , 0 000000 , 2 3 21

0

1 1 21 1 0

1 21 0 21 0 1 2

1 000 , 2

5 000000 , 2 3 0

0 0

0 1

1

= + =

+

−

⋅ +

= ⋅

⋅ + +

=

= + =

⋅ +

−

⋅

⋅ +

= ⋅

⋅

⋅ +

=

= + =

⋅ +

−

⋅

⋅ +

= ⋅

⋅

⋅ +

=

= + =

−

= ⋅

=

+ +

−

⋅ +

⋅

+

⋅ +

−

⋅

⋅ +

⋅

+

⋅ +

−

⋅

⋅ +

−

⋅

h t

k h k z

h z h t f k

h t

k h k z

h z h t f k

h t

k h k z

h z h t f k

t z z t f k

, (30)

{ } { } { } { }

( )

(

0,396450 61⁰^,⁴⁵⁵⁶²¹

)

²^,¹⁹⁵⁹²⁴

31 384615 , 30 333333 1 , 60 500 1 , 0 000000 , 2

1 4 61 1 3 31 1 2 31 1 1 61 0 1

=

⋅ +

⋅

⋅ +

=

⋅ +

⋅

⋅ +

k k

h z

z . (31)

Iteracja druga:

{ }

( )

{ }

(

{ }

) (

^{{ }}

)

( )

{ }

(

{ }

) (

^{{ }}

)

( )

{ }

(

{ }

) (

^{{ }}

)

( ) 0,594280

1 5 , 3

526233 ,

1 0 5 , 3

514135 ,

1 0 5 , 3

453649 ,

1 0 5 , 3

1 500 , 0 500 , 2

5 526233 , 0 500 , 0 195924 , 2 3 1

2 3 2 1

3 1 1 2 4

1 500 , 0 5 , 0 500 , 2

5 514135 , 0 500 , 0 5 , 0 195924 , 2 3 21

1

2 2 21 2 1

2 2 1 1 21 1 2 3

1 500 , 0 5 , 0 500 , 2

5 453649 , 0 500 , 0 5 , 0 195924 , 2 3 21

1

2 1 12 2 1

1 21 1 21 1 2 2

1 500 , 2

5 195924 , 2 3 1

1 1 1 2 1

= + =

+

−

⋅ +

= ⋅

⋅ + +

=

= + =

⋅ +

−

⋅

⋅ +

= ⋅

⋅

⋅ +

=

= + =

⋅ +

−

⋅

⋅ +

= ⋅

⋅

⋅ +

=

= + =

−

= ⋅

=

+ +

−

⋅ +

⋅

+

⋅ +

−

⋅

⋅ +

⋅

+

⋅ +

−

⋅

⋅ +

−

⋅

h t

k h k z

h z h t f k

h t

k h k z

h z h t f k

h t

k h k z

h z h t f k

t z z t f k

, (32)

{ } { } { } { }

( )

(

₆¹⁰^,⁴⁵³⁶⁴⁹ ¹₃⁰^,⁵¹⁴¹³⁵ ¹₃⁰^,⁵²⁶²³³ ¹₆⁰^,⁵⁹⁴²⁸⁰

)

²^,⁴⁵⁶⁶⁴⁶

500 , 0 195924 , 2

1 4 61 1 3 13 1 2 31 1 1 61 0 1

=

⋅ +

⋅

⋅ +

=

⋅ +

⋅

⋅ +

k k

h z

z . (33)

Iteracja trzecia:

{ }

( )

{ }

(

{ }

) (

^{{ }}

)

( )

{ }

(

{ }

) (

^{{ }}

)

( )

{ }

(

{ }

) (

^{{ }}

)

( ) 0,751483

1 5 , 3

674489 ,

1 0 5 , 3

662188 ,

1 0 5 , 3

592484 ,

1 0 5 , 3

1 500 , 0 000 , 3

5 674489 , 0 500 , 0 456646 , 2 3 2

3 3 3 2

3 2 2 3

4

1 500 , 0 5 , 0 000 , 3

5 662188 , 0 500 , 0 5 , 0 456646 , 2 3 12

2

3 2 21 3 2

2 21 2 21 2 3

3

1 500 , 0 5 , 0 000 , 3

5 592484 , 0 500 , 0 5 , 0 456646 , 2 3 21

2

3 2 1 2 1 3

1 12 2 21 2 3

2

1 000 , 3

5 456646 , 2 3 2

2 2 2 3

1

= + =

+

−

⋅ +

= ⋅

⋅ + +

=

= + =

⋅ +

−

⋅

⋅ +

= ⋅

⋅

⋅ +

=

= + =

⋅ +

−

⋅

⋅ +

= ⋅

⋅

⋅ +

=

= + =

−

= ⋅

=

+ +

−

⋅ +

⋅

+

⋅ +

−

⋅

⋅ +

⋅

+

⋅ +

−

⋅

⋅ +

−

⋅

h t

k h k z

h z h t f k

h t

k h k z

h z h t f k

h t

k h k z

h z h t f k

t z z t f k

, (34)

(5)

{ } { } { } { }

( )

(

610,592484 130,662188 130,674489 16⁰^,⁷⁵¹⁴⁸³

)

²^,⁷⁹¹⁴²³

500 , 0 456646 , 2

1 6 4 1 1 3 3 1 1 3 2 1 1 6 1 0 1 1

=

⋅ +

⋅

⋅ +

=

⋅ +

⋅

⋅ +

k k

h z

z . (35)

W tablicach 1 i 2 zestawiono błąd bezwzględny i błąd względny procentowy popełnione przy obli- czaniu wartości funkcji z

( )

t w punkcie t=3,500 przedstawionymi powyżej metodami przy róż- nych długościach kroku ∆t . Wartość rozwiązania analitycznego w tym punkcie, użyta do wyznaczenia błędów bezwzględnego (błąd 1) i względnego procentowego (błąd 2) przedstawionych w ko- lumnach 4 i 5 tablicy, jest równa z^ść

(

3,500

)

=2,791667 .

Tablica 1. Błąd rozwiązania przybliżonego w zależności od długości kroku obliczeń.

Krok Metoda Wartość Błąd 1 Błąd 2 [%] Lof Euler 2,648810 0,142857 5,117265 3 Runge-Kutta II 2,774930 0,016737 0,599534 6 0,5000

Runge-Kutta IV 2,791423 0,000244 0,008740 12 Euler 2,773488 0,018178 0,651160 30 Runge-Kutta II 2,791455 0,000212 0,007594 60 0,0500

Runge-Kutta IV 2,791667 0,000000 0,000000 12000

W ostatniej kolumnie tablicy 1 (Lof) przedstawiono liczbę obliczeń wartości funkcji stojącej po prawej stronie równania (15) koniecznych do wyznaczenia rozwiązania zapisanego w kolumnie trzeciej.

Jak wynika z powyższej tablicy dziesięciokrotne zmniejszenie długości kroku całkowania powodu- je zmniejszenie błędu względnego rozwiązania w metodzie Eulera około dziesięć razy, w metodzie Rungego – Kutty II rzędu około sto razy i w metodzie Rungego – Kutty IV rzędu około dziesięć ty- sięcy razy. Jest to zgodne z stwierdzeniem dotyczącym rzędu każdej z rozważanych metod, zawar- tym na stronie drugiej.

Wyznaczenie wartości rozwiązania przybliżonego w kolejnym punkcie wymaga jednokrotnego obliczenia prawej strony równania (2) w wypadku stosowania metody Eulera, dwukrotnego dla metody Rungego – Kutty II rzędu i czterokrotnego dla metody Rungego – Kutty IV rzędu. W realnych zadaniach obliczanie prawej strony równania (2) jest najbardziej czasochłonną czynnością w każdej iteracji. Wobec tego, aby metoda Rungego – Kutty II rzędu była konkurencyjna w stosunku do metody Eulera, musi dawać wyniki takiej samej dokładności przy dwukrotnie dłuższym kroku obli- czeń, a metoda Rungego – Kutty IV rzędu przy czterokrotnie dłuższym. Jak widać z tablicy 1 (wier- sze 11, 6 i 4) warunki te są spełnione z naddatkiem, gdyż uzyskanie wyniku porównywalnego z wynikiem uzyskanym w trzech krokach metodą Rungego – Kutty IV rzędu (wiersz czwarty tablicy) kosztem dwunastokrotnego obliczenia prawej strony (15) wymaga wykonania około 30 kroków me- todą Rungego – Kutty II rzędu (wiersz szósty tablicy) czyli sześćdziesięciokrotnego obliczenia prawej strony (15) i 3000 kroków metodą Eulera (wiersz jedenasty tablicy) czyli obliczenia prawej strony (15) trzy tysiące razy. Tak więc, spośród trzech prezentowanych metod metoda Rungego – Kutty IV rzędu jest najbardziej efektywna w zastosowaniach praktycznych.

(6)

Rys. 1 Rozwiązanie problemu początkowego (15), (16) w przedziale metodami Eulera (e), Runge- go – Kutty II rzędu (rk2) i Rungego – Kutty IV rzędu (rk4) z krokiem . Dla porównania na wykresie naniesiono również wartości analityczne rozwiązania (an).

W metodach wielokrokowych typu Adamsa – Bashforda i Adamsa – Moultona

2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7

2 3 4 5 6 7 8

e rk2 rk4 an

t