Mariusz PYRZ

SIMR (PW), Instytut Pojazdów

Metody numeryczne w mechanice

Równania różniczkowe zwyczajne

9.

Równania różniczkowe zwyczajne

Równanie różniczkowe pierwszego rzędu

( , )

dy

f x y

dx

=

0

0

(

)

y x

y

=

warunek początkowy

2

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Problem do rozwiązania

: y(x)=?

W równaniach różniczkowych zwyczajnych funkcje niewiadome
zależą od jednej zmiennej niezależnej.

Równania różniczkowe zwyczajne

y

x

y

′ = +

( )

1

x

y x

Ce

x

=

− −

Równanie

Rozwiązanie ogólne

Warunek początkowy

3

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Interpretacja geometryczna

Rysunek 1

(0) 1

y

=

( )

2

1

x

y x

e

x

=

− −

Warunek początkowy

Rozwiązanie szczegόlne

Rząd równania różniczkowego

Dana jest funkcja ciągła f[a,b] x R

n

R

n

f jest klasy C

p

w przedziale

jeżeli jest ciągła i ciągłe są jej pochodne aż do rzędu p

Równanie różniczkowe pierwszego rzędu

x

I

a b

∈ =

0

[ , ]

f

C

I R

p

n

∈

( ,

)

0

′

=

∈ =

∈

∈

y x

f x y x

x

I

a b

y

R

f

C I

( )

( , ( )),

[ , ],

,

( )

0

1

0

4

Równanie różniczkowe rzędu p

Rozwiązanie równania różniczkowego polega na wyznaczeniu funkcji y(x)
spełniających to równanie (jest ich nieskończenie wiele).

′

=

∈ =

∈

∈

y x

f x y x

x

I

a b

y

R

f

C I

( )

( , ( )),

[ , ],

,

( )

0

0

y

x

f x y x y x y

x

y

x

p

p

( )

(

)

( )

( , ( ),

( ),

( ),

,

( ))

=

′

′′

−

…

1

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Warunki początkowe

Wyznaczenie jedynego rozwiązania wymaga podania dodatkowych informacji o

poszukiwanej funkcji – najczęściej są to warunki początkowe.

Niech x

0

oznacza pewien punkt. Warunki początkowe to:

 dla równania

pierwszego

rzędu

wartość funkcji y w punkcie x

0

, czyli

y(x

0

)

 dla równania

rzędu p

5

y(x

0

)

oraz wartości pierwszych (p-1) pochodnych funkcji y w punkcie x

0

czyli

y'(x

0

), y''(x

0

), ... y

(p-1)

( x

0

)

W problemach z zakresu mechaniki opisanych za pomocą równań

różniczkowych warunki początkowe często pojawiają się „naturalnie”.

Uwaga:

Spotyka się również tzw. problemy brzegowe, w których narzucone są

warunki brzegowe y(a) i y(b) na granicach a, b przedziału.

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Zagadnienie Cauchy’ego (zagadnienie początkowe)

Polega na znalezieniu funkcji spełniającej równanie różniczkowe

oraz spełniającej odpowiednio sformułowane warunki początkowe.

Bez warunków początkowych równanie może mieć wiele rozwiązań.
Warunki początkowe pozwalają na odnalezienie jednoznacznego rozwiązania.

y

x

f x y x y x y

x

y

x

p

p

( )

(

)

( )

( , ( ),

( ),

( ),

,

( ))

=

′

′′

−

…

1

6

Styczna do krzywej całkowej y(x) jest określona przez y'(x)=f(x,y(x)).

Każdemu punktowi dziedziny można przyporządkować wektor jednostkowy
o nachyleniu f(x,y(x)) – zbiór takich wektorów określany jest mianem pola.
Scałkowanie równania różniczkowego polega na znalezieniu krzywej y(x)
która ma w każdym punkcie styczną pokrywającą z kierunkiem pola w tym
punkcie.

Interpretacja geometryczna

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Zagadnienie Cauchy’ego (zagadnienie początkowe)

7

Interpretacja geometryczna

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Rysunek 2

Równanie różniczkowe zwyczajne pierwszego rzędu

Zagadnienie Cauchy’ego dla równania różniczkowego pierwszego rzędu
polega na znalezieniu funkcji y(x) klasy C

1

w przedziale I

0

=[a,b]

spełniającej

0

0

0

0

0

( )

( , ( ))

[ , ]

(

)

,

,

y x

f x y x

x

I

a b

y x

y

x y dane

′

=

∈ =

=

8

f jest znaną funkcją,

Tak zdefiniowany problem Cauchy’ego jest równoznaczny z problemem:

Znaleźć y takie, aby spełnić zależność

lub

.

x

I

y

R

0

0

0

∈

∈

,

′

=

z

z

y x dx

f t y t dt

x

x

x

x

( )

( , ( ))

0

0

y

y

f t y t dt

x

x

=

+

z

0

0

( , ( ))

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Istnienie i jednoznaczność rozwiązania

Jeżeli funkcja f jest ciągła i spełnia warunek Lipschitza dla y, to istnieje
funkcja y klasy C

1

w I

0

która stanowi jedyne rozwiązanie problemu

Cauchy’ego.

Przypomnienie:

Funkcja f(x,y) spełnia warunek Lipschitza dla y w I

0

x R jeżeli istnieje

stała

taka, że

0

< ∈

K

R

9

Oznacza to, że szybkość zmian funkcji jest ograniczona, Moglibyśmy
zamiast spełnienia warunku Lipschitza zażądać aby f była klasy C

1

, ale

taki warunek jest zbyt surowy gdyż istnieją funkcje f które mają nieciągłe
pochodne ale spełniają warunek Lipschitza.

0

< ∈

K

R

f x y

f x y

K y

y

x

a b

y y

R

n

( ,

)

( ,

)

[ , ],

,

1

2

1

2

1

2

−

<

−

∀ ∈

∀

∈

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Rozwiązywanie numeryczne – dyskretyzacja

Przybliżone wartości poszukiwanej funkcji y(x

i

) obliczane są

w n punktach x

i

z przedziału I

0

=[a,b].

Zasada:

Dzielimy przedział I

0

za pomocą n+1 punktów x

0

,x

1

, ... x

n

(zazwyczaj równoodległych)

x

0

=a, x

i+1

=x

i

+h, 0<=i<=n gdzie krok dyskretyzacji h=(b-a)/n

Obliczamy n wartości y

1

, y

2

, …, y

n

w punktach x

i

tj. przybliżone wartości y(x

1

), y(x

2

),

10

Obliczamy n wartości y

1

, y

2

, …, y

n

w punktach x

i

tj. przybliżone wartości y(x

1

), y(x

2

),

…, y(x

n

) za pomocą ogólnego wzoru y

i+1

~y

i

+F(x

i

, y

i

,h)

Łączymy otrzymane punkty aby uzyskać przebieg funkcji y(x) w I

0

(

można wykorzystać procedury interpolacji lub aproksymacji)

Szacujemy bład dyskretyzacji (który zależy od h)

e

i

=y(x

i

)-y

i

(y(x

i

) – wartość dokładna, y

i

– wartość przybliżona)

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Rozwiązanie numeryczne

11

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Rysunek 3

Rodzaje metod

Metody o pojedynczych krokach (jednokrokowe)

y

i+1

jest obliczana w funkcji wartości y

i

, x

i

i pozostałych danych

Metody o pojedynczych krokach (explicite) konstruują ciąg przybliżeń

posługując się zależnością

(funkcja F może być nieliniowa względem f(x) )

y

y x

i

N

i

i

≈

=

( ),

, ,

,

0 1 …

y

y x

y

y

hF x y h

i

N

i

i

i

i

0

0

1

0 1

=

= +

=

+

(

)

( ,

, ),

, ,

,

…

12

(funkcja F może być nieliniowa względem f(x) )

Metody o połączonych krokach (wielokrokowe)

y

i+1

jest obliczana na podstawie wartości y

i

,y

i-1

,y

i-2

,… obliczonych wcześniej.

Metody liniowe można zapisać:

Metoda nazywamy k – krokową jeżeli

α

α

α

β

β

β

k

i k

i

i

k

i k

i

i

j

j

j

j

j

y

y

y

h

f

f

f

i

N

k

y

y h

j

k

f

f x y

+

+

+

+

+ +

+

=

+ +

+

=

−

=

=

−

=

…

…

…

…

1

1

0

1

1

0

0 1

0 1

1

(

),

, ,

,

( ),

, ,

,

,

(

,

)

α

α

β

k

et

≠

+

≠

0

0

0

0

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Stabilność

Pojęcie stabilności metody oznacza, że małe zaburzenie danych y

0

(warunków początkowych) i funkcji F pociąga za sobą jedynie małe
zaburzenie rozwiązania i to niezależnie od wielkości kroku h.

Zbieżność

Pojęcie zbieżności oznacza, że przybliżone rozwiązanie powinno dążyć do
rozwiązania dokładnego w każdym punkcie przedziału I

0

wraz ze

13

rozwiązania dokładnego w każdym punkcie przedziału I

0

wraz ze

zmniejszaniem kroku dyskretyzacji h (tj. h dążącym do zera).

Metoda zdefiniowana przez schemat

jest zbieżna jeśli dla każdego rozwiązania dokładnego y(x) ciąg y

i

spełnia

warunek

dla każdego warunku początkowego.

wyraża błąd globalny ciągu y

i

w stosunku do rozwiązania

dokładnego y(x

i

) (ważne w praktyce)

0

0

1

(

)

( ,

, ),

0,1,

,

i

i

i

i

y

y x

y

y

hF x y h

i

N

+

=

= +

=

…

lim sup

( )

h

i

i

y

y x

→

−

=

0

0

sup

( )

y

y x

i

i

−

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Spójność schematu rozwiązania

Pojęcie spójności określa w jaki sposób schemat rozwiązanie odzwierciedla
równanie różniczkowe y'=f(x,y).

Metoda zdefiniowana za pomocą schematu rozwiązania

jest określana jako spójna z równaniem różniczkowym y'=f(x,y), y

0

=y(x

0

) jeżeli

dla każdego rozwiązania y problemu zachodzi

lub

lim sup

(

)

( )

( , ( ))

i

i

y x

y x

F x y x

+

−

−

L

M

O

P

=

1

0

lim

(

)

0

y x

y

ε

=

−

=

∑

∑

0

0

1

(

)

( ,

, ),

0,1,

,

i

i

i

i

y

y x

y

y

hF x y h

i

N

+

=

= +

=

…

14

lub

Twierdzenie: Schemat jest spójny wtedy i tylko wtedy jeżeli

Jeżeli funkcja F spełnia warunki F(x,y,0)=f(x,y) i jeżeli F spełnia warunek
Lipschitza, to metoda jest spójna.

Jeżeli metoda jest spójną i stabilna to jest zbieżna

.

lim sup

(

)

( )

( , ( ))

h

i n

i

i

i

i

y x

y x

h

F x y x

→

≤ ≤ −

+

−

−

L

NM

O

QP

=

0 0

1

1

0

1

1

0

0

1

0

1

lim

(

)

0

i

i

i

h

i n

i n

y x

y

ε

+

+

→

≤ ≤ −

≤ ≤ −

=

−

=

∑

∑

F x y

f x y

x

I

y

R

( , , )

( , )

0

0

=

∀ ∈

∀ ∈

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Rząd metody

Rząd metody wskazuje ilościowo w jaki sposób metoda jest zbieżna w
zależności od h, czyli podaje informacje o sposobie w jaki

zmierza do 0 wraz z h.

Metoda jest rzędu p (p > 0) jeżeli dla każdego rozwiązania y(x) istnieje K<>0

takie że:

ou

sup

( )

y

y x

i

i

−

sup

0

1

≤ ≤ −

≤

i n

n

p

Kh

ε

sup

(

)

( )

( , ( ), )

0

1

1

≤ ≤ −

+

−

−

≤

i n

i

i

i

i

p

y x

y x

h

F x y x

h

Kh

15

gdzie K jest niezależna od h, ale zależy od funkcji y i F, o których zakładamy
że są różniczkowalne w wystarczającym stopniu.

ε

n

(x) jest zatem rzędu p ze względu na h (

ε

n

(x)

∈

O(h

p

)).

Dany schemat rzędu p jest oczywiście spójny.

Metoda zbiega sie tym szybciej, im wyższy jest jej rząd p.

Jeżeli na przykład zmniejszymy dwukrotnie krok, to nakład pracy wzrasta dwukrotnie ale
oszacowanie błędu zmniejsza się 2p – krotnie. Błąd lokalnie zmniejsza sie jak h

p

.

Błąd metody Euler’a-Cauchy’ego jest proporcjonalny do h.

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Metoda Eulera (stycznej)

Krzywa przechodząca przez (x

i

, y

i

) jest zastępowana przez jej styczną

Rozwiniecie Taylor funkcji y w otoczeniu punktu x

i

R

T

– reszta

z rozwinięcia

R

T

/h wystarczająco mały y’(x

i

) ~ [y(x

i+1

)- y(x

i

)]/h

Punkt o współrzędnych (x , y ) jest położony na stycznej w punkcie (x , y )

y x

y x

hy x

R

i

i

i

T

(

)

( )

( )

+

=

+ ′

+

1

′

=

−

−

+

y x

y x

y x

h

R

h

i

i

i

T

( )

(

)

( )

1

y x

y x

hf x y x

i

i

i

i

(

)

( )

( , ( ))

+

=

+

1

16

Punkt o współrzędnych (x

i+1

, y

i+1

) jest położony na stycznej w punkcie (x

i

, y

i

)

do krzywej całkowej przechodzącej przez ten punkt

Schemat (algorytm):

Funkcja F(x,y,h) charakteryzuje ogólnie schemat rozwiązania. Metoda Euler’a-Cauchy’ego jest
przypadkiem szczególnym z funkcją F(x,y,h)=f(x,y(x)) niezależną od h.

Metoda ta jest zazwyczaj wykorzystywana aby oszacować przebieg funkcji – następnie
można przejść do metody wyższego rzędu.

0

0

1

(

),

( ,

, )

( ,

, )

( , ( ))

i

i

i

i

i

i

i

i

y

y x

y

y

hF x y h

F x y h

f x y x

+

=

= +

=

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Metoda Eulera

17

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Rysunek 4

Zmodyfikowana metoda Eulera

Schemat :

Wykorzystujemy punkt pośredni

zdefiniowany jako

położony na stycznej w punkcie (x

i

,y

i

) do krzywej całkowej przechodzącej

0

0

1

( )

( ,

, )

( ,

, )

(

,

( ,

))

2

2

i

i

i

i

i

i

i

i

i

i

y

y x

y

y

hF x y h

h

h

F x y h

f x

y

f x y

+

=

= +

=

+

+

(

,

)

/

/

x

y

i

i

+

+

1 2

1 2

x

x

h

y

h

f x y

i

i

i

i

i

+

+

= +

=

1 2

1 2

2

2

/

/

,

( ,

)

18

i

i

przez ten punkt

Otrzymujemy następujący algorytm obliczania y

i+1

Prosta przechodząca przez punkty (x

i

,y

i

) , (x

i+1

, y

i+1

) jest równoległa do stycznej

w punkcie

do krzywej całkowej przechodzącej przez ten punkt

y x

y

x

x

h

y

h

f x y

y

y

hf x

y

i

i

i

i

i

i

i

i

i

(

)

,

,

( ,

),

(

,

)

/

/

/

/

0

0

1 2

1 2

1

1 2

1 2

2

2

=

= +

=

= +

+

+

+

+

+

(

,

)

/

/

x

y

i

i

+

+

1 2

1 2

M.Pyrz Metody numeryczne w mechanice – Równania różniczkwe zwyczajne 04.2012

Zmodyfikowana metoda Eulera

19

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Rysunek 5

Zmodyfikowana metoda Eulera (2)

Schemat:

Wykorzystujemy punkt pośredni (x

i+1

, y

i+1/2

) zdefiniowany jako

położony dla odciętej x na stycznej w punkcie (x ,y ) do krzywej całkowej

0

0

1

( )

( ,

, )

1

( ,

, )

[ ( ,

)

(

,

( ,

))]

2

i

i

i

i

i

i

i

i

i

i

i

i

y

y x

y

y

hF x y h

F x y h

f x y

f x

h y

hf x y

+

=

= +

=

+

+

+

y

y

hf x y

i

i

i

i

+

= +

1 2

/

( ,

)

20

położony dla odciętej x

i+1

na stycznej w punkcie (x

i

,y

i

) do krzywej całkowej

przechodzącej przez ten punkt

Uzyskujemy następujący algorytm obliczania y

i+1

Prosta przechodząca przez punkty (x

i

,y

i

) , (x

i+1

, y

i+1

) jest równoległa do

stycznej do krzywej całkowej przechodzącej przez punkt (x

i+1/2

, y

i+1/2

).

0

0

1/ 2

1

1

1/ 2

( )

( ,

)

[ ( ,

)

(

,

)]

2

i

i

i

i

i

i

i

i

i

i

y x

y

y

y

hf x y

h

y

y

f x y

f x

y

+

+

+

+

=

= +

= +

+

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Zmodyfikowana metoda Eulera (2)

21

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Rysunek 6

Metoda Heun’a

Schemat:

Algorytm obliczeń:

y

y x

y

y

hF x y h

F x y h

f x y

f x

h y

hf x y

i

i

i

i

i

i

i

i

i

i

i

i

0

0

1

1

4

3

4

2

3

2

3

=

= +

=

+

+

+

+

(

)

( ,

, )

( ,

, )

[ ( ,

)

(

,

( ,

))]

22

y

y x

x

x

h

y

y

hf x y

y

y

h

f x y

f x

y

i

i

i

i

i

i

i

i

i

i

i

i

0

0

2 3

2 3

1

2 3

2 3

2

3

2

3

4

3

4

=

= +

= +

= +

+

+

+

+

+

+

(

)

( ,

)

[ ( ,

)

(

,

)]

/

/

/

/

Metody Rungego - Kutty

Metody Runge-Kutta (zwane również RK

R

, gdzie R jest liczbą całkowitą

określającą rząd metody) stanowią klasę metod jednokrokowych opartych
na następującej zasadzie:
- W każdej iteracji, w punktach innych niż te dla których poszukujemy
rozwiązania, obliczane są wartości pośrednie (x

nj

, y

nj

) wykorzystując

y

y

y

y

h

f x

y

j

R

n

n

nj

n

jk

nk

nk

j

0

1

1

=

=

+

=

−

∑

α

(

,

),

,

,

…

23

- jako rozwiązanie przyjmujemy y

n+1

= y

nR

.

Wartości 

jk

, 

k

są dobierane w taki sposób, aby rozwinięcie y

nj

dla

kolejnych rosnących potęg h przybliżało najlepiej jak można rozwinięcie w
szereg Taylor rozwiązania problemu.

y

y

h

f x

y

j

R

x

x

h

k

R

nj

n

jk

nk

nk

k

nk

n

k

k

R

0

0

1

0 1

1

0

1

=

+

=

=

+

∈

=

=

=

=

∑

α

θ

θ

θ

θ

(

,

),

,

,

,

[ , ],

,

,

…

…

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Metody Rungego - Kutty

Metodę Euler’a-Cauchy’ego można traktować jako metodę RK

1

gdzie R=1,



10

=1 

0

=0 :

Metoda stycznej ulepszonej: 

10

=1, 

20

=0, 

21

=1, 

0

=0, 

1

=1/2, 

2

=1.

Wartości pośrednie są zdefiniowane jako

y

y

y

h

f x y

y

h f x y

n

n

n

n

n

n

n

n

+

=

=

+

=

+

1

1

10

0

α

(

,

)

(

,

)

y

y

y

y

h

f x y

y

y

y

hf x

h

y

n

n

n

n

n

n

n

n

n

n

n

0

1

1

2

1

2

2

=

=

+

=

=

+

+

+

,

(

,

),

(

,

)

24

Metoda Euler zmodyfikowana: 

10

=1, 

20

=1/2, 

21

=1/2, 

0

=0, 

1

=1, 

2

=1.

Wartości pośrednie są zdefiniowane za pomocą zależności

Metoda Heuna: 

10

=2/3, 

20

=1/4, 

21

=3/4, 

0

=0, 

1

=1, 

2

=1.

n

n

n

n

n

n

n

n

n

n

n

0

1

1

2

1

2

2

+

y

y

y

y

h

f x

y

y

y

y

h

f x

y

f x

h

y

n

n

n

n

n

n

n

n

n

n

n

n

n

0

1

1

2

1

2

2

2

=

=

+

=

=

+

+

+

+

,

(

,

),

[ (

,

)

(

,

)]

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Metoda Rungego – Kutty czwartego rzędu (klasyczna)

Schemat:

gdzie

stosunkowo kosztowna gdyż

y

y

hF x y h

F x y h

k

k

k

k

n

n

n

n

n

n

+

=

+

=

+

+

+

1

1

2

3

4

1

6

2

2

(

,

, )

(

,

, )

(

)

k

f x y

k

f x

h

y

h

k

n

n

1

=

=

+

+

(

,

)

(

,

)

25

stosunkowo kosztowna gdyż

potrzebuje wielu obliczeń funkcji

f(x,y) w każdym kroku

Metoda jest 4-go rzędu i odpowiadają jej parametry



10

=1/2, 

20

=0, 

21

=1/4, 

30

=0, 

31

=0, 

32

=1, 

40

=1/6,



41

=1/3, 

42

=1/3, 

43

=1/6, 

0

=0, 

1

=1/2, 

2

=1/2, 

3

=1

k

f x

h

y

h

k

k

f x

h

y

h

k

k

f x

h y

hk

n

n

n

n

n

n

2

1

3

2

4

3

2

2

2

2

=

+

+

=

+

+

=

+

+

(

,

)

(

,

)

(

,

)

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Metody wielokrokowe

Pozwalają wyeliminować niedogodność liczenia wartości f(x,y) w punktach
pośrednich, nie służących bezpośrednio jako rozwiązanie

Schemat :

nie występują wartości pośrednie

(wykorzystuje się wartości obliczone na poprzednich etapach)

Metody te nie są samo startujące

y

F x y h

x

y

h

n

n

n

n

n k

n k

n k

+

−

−

−

=

1

(

,

,

,

,

,

,

)

…

26

Metody te nie są samo startujące

aby zacząć iteracje niezbędne jest zastosowanie metody explicite

umożliwiającej obliczenie pierwszych wartości

Wzór ogólny:

y

n+i

są dane lub obliczone wcześniej, k jest liczbą kroków

wykorzystywanych do obliczeń, 

i

,

i

– stałe niezależne od n

0

0

(

,

),

0

,

0

k

k

i

n i

i

n i

n i

k

i

i

y

h

f x

y

n

N

k

α

β

α

+

+

+

=

=

=

≤ ≤ −

≠

∑

∑

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Metody wielokrokowe

27

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Rysunek 7

Metody wielokrokowe

Jeżeli 

k

=0 to metoda jest jawna („explicite”)

(można bezpośrednio otrzymać y

n+k

jako funkcję y

n

,y

n+1

,..., y

n+k-1

)

Jeżeli 

k

0 to metoda jest niejawna („implicite”)

(y

n+k

jest zdefiniowane za pomocą równania niejawnego w postaci

y

n+k

=f(x

n+k

, y(x

n+k

)) + „znane człony” )

Rząd metod wielokrokowych

28

Rząd metod wielokrokowych

Metoda k-krokowa jest rzędu p jeżeli, niezależnie od rozwiązania y, mamy

gdzie C nie zależy od h.

Przykład:

metoda Adams’a-Basforth’a, Adams’a-Moulton’a,
metoda predykcji-korekcji Adamsa rzędu 4

sup

(

)

(

, (

))

0

0

0

1

≤ ≤ −

+

+

+

=

=

−

≤

∑

∑

n N k

i

n i

i

n i

n i

i

k

i

k

p

h

y x

f x

y x

Ch

α

β

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Metody prognozy i korekcji

Idea polega najpierw na zastosowaniu sformułowania jawnego pozwalającego
na oszacowanie pierwszego przybliżenia wartości y

i

(

prognoza

)

Następnie stosuje się sformułowanie niejawne aby „poprawić" wartość y

i

(

korekcja

), wykorzystaną w drugim członie poprzedniego wyrażenia

Metoda typu "prognoza-korekcja" drugiego rzędu

1

1

1

( ,

,...,

,

)

i

i

i

i

i

i

y

y

F x x

y y

+

+

+

= +

*

*

1

( ,...,

)

i

i

i

i

y

y

F x

y

+

= +

29

Metoda typu "prognoza-korekcja" drugiego rzędu

1. Wybrać krok h (y

0

jest dane).

2. Obliczyć wartość y

1

stosując np. metodę RK

2

3. Obliczać kolejne wartości y

2

, y

3

,…(n=1,2,…) stosując wzory

*

1

1

1

*

1

1

1

3

(

,

)

(

,

)

2

2

(

,

)

(

,

)

2

n

n

n

n

n

n

n

n

n

n

n

n

h

h

y

y

f x y

f x

y

h

y

y

f x y

f x

y

+

−

−

+

+

+

=

+

−





=

+

+





1

n

n

x

x

h

+

= +

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Metoda typu "prognoza-korekcja" Adams’a rzędu 4

1. Wybrać krok h (y

0

dane).

2. Obliczyć 3 pierwsze wartości y

1

, y

2

, y

3

stosując metodę samo startującą –

np. RK

4

(n=0,1,2) :

k

f x y

k

f x

h

y

h

k

k

f x

h

y

h

k

k

f x

h y

hk

n

n

n

n

1

2

1

2

2

=

=

+

+

=

+

+

=

+

+

(

,

),

(

,

),

(

,

),

(

,

)

y

y

h

k

k

k

k

x

x

h

n

n

n

n

+

+

=

+

+

+

+

=

+

1

1

2

3

4

1

6

2

2

(

),

30

3. Obliczać następne wartości y

4

, y

5

,…(n=3,4,…) stosując wzory

k

f x

h

y

h

k

k

f x

h y

hk

n

n

n

n

3

2

4

3

2

2

=

+

+

=

+

+

(

,

),

(

,

)

y

y

h

f x y

f x

y

f x

y

f x

y

x

x

h

y

y

h

f x

y

f x y

f x

y

f x

y

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

n

+

−

−

−

−

−

−

+

+

+

+

−

−

−

−

=

+

−

+

−

=

+

=

+

+

−

+

1

1

1

2

2

3

3

1

1

1

1

1

1

2

2

24

55

59

37

9

24

9

19

5

*

*

(

,

)

(

,

)

(

,

)

(

,

) ,

,

(

,

)

(

,

)

(

,

)

(

,

)

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Układy równań różniczkowych pierwszego rzędu

Układu p równań różniczkowych pierwszego rzędu

po wprowadzeniu wektorów

′

=

′

=

′

=

y

x

f

x y

x y

x

y

x

y

x

f

x y

x y

x

y

x

y

x

f

x y

x y

x

y

x

p

p

p

p

p

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( ,

( ),

( ),

,

( ))

( )

( ,

( ),

( ),

,

( ))

( )

( ,

( ),

( ),

,

( ))

1

1

1

2

2

2

1

2

1

2

…

…

⋯

…

L O

′

L O

L

O

y

y

f

x y

x y

x

y

x

( ,

( ),

( ),

,

( ))

…

31

można zapisać

warunki początkowe

y

y

f

y

=

L

N

M

M

M

M

O

Q

P

P

P

P

′ =

′

′

′

L

N

M

M

M

M

O

Q

P

P

P

P

=

L

N

M

M

M

M

O

Q

P

P

P

P

y

y

y

y

y

y

x

f

x y

x y

x

y

x

f

x y

x y

x

y

x

f

x y

x y

x

y

x

p

p

p

p

p

p

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( , )

( ,

( ),

( ),

,

( ))

( ,

( ),

( ),

,

( ))

( ,

( ),

( ),

,

( ))

1

2

1

2

1

1

2

2

1

2

1

2

⋮

⋮

…

…

⋮

…

y

f

y

f

y

( )

,

:( , )

( , )

(

)

x

R

x

x

R

R

R

p

p

p

∈

→

×

→

′

=

=

∈

y

f

y

y

y

y

0

0

( )

( , ( ))

(

)

,

,

x

x

x

x

donné

x

I

0

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Rozwiązanie układu równań różniczkowych pierwszego rzędu

Całkowanie układu równań różniczkowych można przeprowadzić podobnie jak
rozwiązuje się pojedyncze równanie pierwszego rzędu, stosując
przedstawione wcześniej metody

Metoda Eulera:

gdzie dla i=1,2,…,p :

(

)

(

)

(

, (

) , (

) ,

, (

) )

( )

( )

( )

( )

( )

( )

y

y

hf

x

y

y

y

i

n

i

n

i

n

n

n

p

n

+

=

+

1

1

2

…

y

y

f

y

n

n

n

n

h

x

+

=

+

1

(

,

)

32

Zmodyfikowana metoda Eulera:

gdzie dla i=1,2,…,p :

-

wartości i-tej składowej y

n

i f=y'

n

.

y

y

f

y

n

n

n

n

h

x

+

=

+

1

(

,

)

y

y

f

y

f

y

n

n

n

n

n

n

h

x

x

+

+

+

=

+

+

1

1

1

2

(

,

)

(

,

)

( )

1

( )

( )

(1)

( )

( )

1

(1)

1

( )

1

(

)

(

)

[

(

, (

) ,..., (

) )

(

, (

)

,..., (

)

)]

2

i

n

i

n

i

n

n

p

n

i

n

n

p

n

h

y

y

f

x

y

y

f

x

y

y

+

+

+

+

=

+

+

(

) ,

(

, (

) ,

, (

) )

( )

( )

( )

( )

y

f

x

y

y

i

n

i

n

n

p

n

1

…

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

• Metoda RK

4

:

y

y

k

k

k

k

n

n

n

n

n

n

+

=

+

+

+

+

1

1

2

3

4

1

6

2

2

(

)

(

)

(

)

(

)

(

)

(

,

)

(

)

(

,

(

) )

(

)

(

,

(

) )

k

f

y

k

f

y

k

k

f

y

k

1

2

1

3

2

2

1

2

2

1

2

n

n

n

n

n

n

n

n

n

n

n

h

x

h

x

h

h

x

h

=

=

+

+

=

+

+

Rozwiązanie układu równań różniczkowych pierwszego rzędu

33

(

)

(

,

(

) )

k

f

y

k

4

3

2

2

1

2

n

n

n

n

h

x

h

=

+

+

(

)

(

,

)

(

,

)

(

,

)

( )

( )

( )

k

y

y

y

1

1

2

n

n

n

n

n

p

n

n

h

f

x

f

x

f

x

=

L

N

M

M

M

M

O

Q

P

P

P

P

⋯

(

)

(

,

(

) )

(

,

(

) )

(

,

(

) )

( )

( )

( )

k

y

k

y

k

y

k

2

1

1

2

1

1

2

1

2

2

1

2

2

1

2

n

n

n

n

n

n

n

p

n

n

n

h

f

x

h

f

x

h

f

x

h

=

+

+

+

+

+

+

L

N

M

M

M

M

M

M

M

O

Q

P

P

P

P

P

P

P

⋯

⋯

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Równanie różniczkowe wyższego rzędu

Za pomocą zmiany zmiennych przekształcamy równanie różniczkowe rzędu p > 1
do układu p równań różniczkowych rzędu 1.

Niech będzie dane równanie różniczkowe rzędu p (p>1)

Wprowadzamy nowe zmienne

y

f x y x y x

y

x

p

p

( )

(

)

( , ( ),

( ),

,

( ))

=

′

−

…

1

y x

y x

( )

( )

=

′ =

y

y

1

2

Do rozwiązania:
Układ p równań
różniczkowych 1-go

34

y x

y x

y x

y x

y

y x

y

x

y

y

x

y

x

y

p

p

p

1

2

1

3

2

1

1

( )

( )

( )

( )

( )

( )

( )

( )

(

)

=

= ′

= ′

= ′′

= ′

=

= ′

−

−

⋯

′ =
′ =

′ =
′ =

−

y

y

y

y

y

y

y

f

p

p

p

1

2

2

3

1

⋯

różniczkowych 1-go
rzędu w którym
niewiadomymi są
funkcje

y

1

, y

2

, …, y

p

.

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012

Równanie różniczkowe wyższego rzędu
Przekształcenie warunków początkowych

Warunki początkowe

Nowe warunki początkowe

0

0

0

0

( )

( )

y x

y

y x

y

=

′

′

=

1

0

0

2

0

0

( )

(

)

y x

y

y x

y

=

′

=

35

0

0

(

1)

(

1)

0

0

( )

(

)

p

p

y x

y

y

x

y

−

−

′′

′′

=

=

…

3

0

0

(

1)

0

0

( )

( )

p

p

y x

y

y

x

y

−

′′

=

=

…

M.Pyrz Metody numeryczne w mechanice – Równania różniczkowe zwyczajne 04.2012