METODY NUMERYCZNE

METODY NUMERYCZNE

dr hab.inż. Katarzyna Zakrzewska, prof.AGH, Katedra Elektroniki, AGH

e-mail: zak@agh.edu.pl

Wykład 7

Równania różniczkowe - przegląd

Równania różniczkowe - wprowadzenie

Równania różniczkowe są popularnie spotykane we

wszystkich dziedzinach nauk ścisłych i przyrodniczych a szczególnie w:

• Fizyce (np. równania Maxwell’a)

• Mechanice (np. równania ruchu harmonicznego)

• Elektronice (np. stany nieustalone w obwodach elektrycznych)

• Automatyce (np. warunki sterowalności układu)

• i wielu innych dziedzinach nauki i techniki

Równania różniczkowe zwyczajne

Równania różniczkowe zwyczajne – jest to równanie w którym występują stałe oraz funkcje niewiadome i pochodne funkcji niewiadomych zależne od jednej zmiennej niezależnej.

Przykład:

) (

... _{2 2 1}

2 1 3

1

x F y

dx k k dy dx

y k d

dx y k d

dx y d

n n n n

n  ^_     

K b

x M

2

0

2

  Kx 

dt b dx dt

x

M d

Cząstkowe równania różniczkowe

Cząstkowe równanie różniczkowe – jest to równanie zawierające funkcję niewiadomą dwóch lub więcej zmiennych oraz

niektóre z jej pochodnych cząstkowych.

Jednym z najprostszych równań różniczkowych cząstkowych

jest równanie transportu:

 

 











 













n t

x u x

u u

x t u u

u b

u

,..., )

, (

0

1

Cząstkowe równania różniczkowe

0 )

, (

) ,

) ( (

) ,

( )

(























sb x

s t

u sb

x s t

ds u s dz

sb x

s t

u s

z

x t

Zauważamy, że pochodna kierunkowa funkcji u w kierunku wektora v=(1,b) є Rⁿ⁺¹ znika. Zatem ustalając dowolny punkt

(t,x) є R₊ x Rⁿ i kładąc dla s є R dostajemy:

Zatem z(s) jest funkcją stałą.

Ustalając wartość rozwiązania na każdej prostej równoległej do wektora (1,b) dostajemy rozwiązanie zadania.

Cząstkowe równania różniczkowe – zagadnienia początkowe



 



















n

n

R x

dla x

g x

u

R x

t dla u

b u

) ( )

, 0 (

, 0 0

Zagadnienia początkowe, zakładamy, że w chwili t=0 zadana jest wartość funkcji u(0,x). Wówczas zagadnienie

początkowe:

ma rozwiązanie:

R

x t

dla tb

x g x

t

u ( , )  (  )  0 , 

Jeśli funkcja g jest klasy C¹ to rozwiązanie równania jest rozwiązaniem klasycznym oraz jest ono jednoznaczne

Rozwiązywanie zagadnień początkowych

Wprowadzimy teraz kilka oznaczeń, niech:

 Y(x) – oznacza dokładne rozwiązanie

 y(x) – oznacza rozwiązanie przybliżone

N i

b x

x

gdzie

f y

x f

y x

y y

Y x

dx f x Y dY

x Y

Y

i i

i i

i i

i i i

i i

i

,...., 2

, 1 ,

; (

:

) ,

( ),

(

) ,

) ( ), (

(

' '

















są punktami, w których wyznaczamy przybliżone

Błąd metody

Wielkość Tn nazywamy błędem metody powstałym przy przejściu od x_n do x_n+1

ih x

x



 h – krok całkowania

Błąd metody możemy wyrazić jako funkcję zmiennej h i przedstawić w postaci:

) (

1 







n

h O h

T 

Gdzie stała , to liczbę p będziemy nazywać rzędem metody przybliżonej

 0



Cząstkowe równania różniczkowe – zagadnienia niejednorodne



 



















n

n

R x

dla x

g x

u

R x

t dla f

u b

u

) ( )

, 0 (

, 0

W celu rozwiązania zagadnienia niejednorodnego:

podstawmy:

wówczas:

) ,

(

) ,

( )

, ) (

(

bs x

s t

f

b bs

x s t

u bs

x s t

ds u s dz

t

























) ,

( )

(s u t s x sb

z   

Cząstkowe równania różniczkowe – zagadnienia niejednorodne







































t t

t

ds b

t s x

s f

ds sb x

s t

f

ds t dz

z z

tb x

g x

t y

0 0

0

) ) (

, (

) ,

(

) ( )

0 ( )

( )

,

Zatem:

(

Rozwiązaniem zagadnienia jest więc:







 g x tb ^t f s x s t b ds x

t

u( , ) ( ) 0 ( , ( ) )

Całka zupełna dla równań rzędu 1

0 )

, ,

( x u u  F

Ogólne równanie różniczkowe cząstkowe rzędu 1 można zapisać w postaci:

gdzie:

 

R R

R F

u u

u

R u

x

n x x























 :

,..., :

,

2 1

Całka zupełna dla równań rzędu 1 – zagadnienia brzegowe



















i powierzchn na

g u

R obszarze

w u

u x

F ( , , ) 0

Ogólne równanie różniczkowe cząstkowe rzędu 1 można

zapisać w postaci, spróbujemy rozwiązać warunki brzegowe (zagadnienie Cauchy’ego)

Zakładamy, że funkcje F i g oraz powierzchnia Г są dostatecznie gładkie.

Całka zupełna dla równań rzędu 1 – zagadnienia brzegowe



 

 







j i j

ij j

ij i

x x

u u n i

dla s

x s

x u

s p

1

2

,..., 1

) ( ))

( ( )

(

Staramy się połączyć punkt xєΩ z pewnym punktem x₀єГ pewną krzywą ɣ w taki sposób, aby można było policzyć wartości

rozwiązania u wzdłuż tej krzywej:

Chcemy dobrać krzywą x(s) tak, aby można było policzyć z(s) i p(s). W tym celu policzymy dp(s)/d(s):

 

)) (

( )

( )),

( ( )

(

), (

s x u s

p s

x u s

z

R I

s s

x











 

Całka zupełna dla równań rzędu 1 – zagadnienia brzegowe

n j

dla s

p s

z s

p x s F

x

j

j

( ) ( ( ), ( ), ( ))  1 ,...,



 

Z drugiej strony różniczkując równanie ogólne względem x_i:

Zakładamy, że:







 

 

 

 

 



j

ij j

x i

n i

u u u

p x u F

u u

z x u F

u x x

F

i 1

,...

1 0

) ,

, ( )

, , ( )

,

,

(

Całka zupełna dla równań rzędu 1 – zagadnienia brzegowe

 

 





 







j

n

j j

j j

j

s p s

z s

p x s F

p s

x s

x x s u

z

1 1

)) ( ), ( ), ( ( )

( )

( ))

( ( )

(

Otrzymujemy wtedy:

Ostatecznie otrzymujemy:





z s F

p s

z s x x

s F

p ⁱ

i

i( ) ( ( ), ( ), ( )) ( ), ( ), ( ) ( ) 1,...,



 



 

Całka zupełna dla równań rzędu 1 – zagadnienia brzegowe

)) ( ), ( ), ( ( )

(

) ( ))

( ), ( ), ( ( )

(

) ( ))

( ), ( ), ( ( ))

( ), ( ), ( ( )

(

s p s

z s

x F s

x

s p s

p s

z s x F s

z

s p s

p s

z s

x F s

p s

z s x F s

p

p p

z x





















Stosując zapis wektorowy otrzymujemy układ (2n+1) równań różniczkowych zwyczajnych zwany układem charakterystyk całki zupełnej.

Całka zupełna dla równań rzędu 1 – zagadnienia brzegowe - przykład

 

 

 

























0 ,

0 )

( )

, (

0 ,

0

2 1

1 2

1

2 1

1 2 2

1

x x

x g x

x u

x x

u u

x u

x

Rozpatrzmy układ:

Wówczas równania charakterystyk mają postać:

z z

x x

x

x

 

,



,



Całka zupełna dla równań rzędu 1 – zagadnienia brzegowe - przykład

0 2 ,

0

) (

) ( , sin )

( ,

cos )

(

0

0 0

2 0

1

 











s x

e x

g s

z s x

s x

s x

s

x

Rozwiązując ten układ równań z uwzględnieniem warunku brzegowego dostajemy

Ostatecznie rugując parametr s dostajemy rozwiązanie:

0 ,

1 ,

exp )

( )

,

(

1 2 2

2 2

1 2

1

 

 



 

 

 x x

x arctg x x

x g

x

x

u

Równanie liniowe rzędu 2

 

 





n



l k

n

k

xk k

xl xk l

k

x u x b x u x c x u f x

a

1

, 1

,

,

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

Ogólne równanie różniczkowe cząstkowe liniowe rzędu 2 określone w obszarze Ω ↄ Rⁿ ma postać:

Ogólne równanie cząstkowe drugiego rzędu dwóch zmiennych niezależnych liniowe możemy zapisać:

) ,

, , , ( )

, ( )

, ( 2 )

,

(

2 2

2 2

y f x

f f y x y F

y f x xy C

y f x x B

y f x

A 





 



 



 





A,B,C są określone w pewnym obszarze Ω ↄ R² a F jest

Transformacja Laplace’a

) ( )

(

0

s F dt

t f

e









Przy czym o funkcji f(t) zakładamy że jest funkcją rzeczywistą zmiennej rzeczywistej t określoną dla każdej wartości t > 0 i przedziałami ciągłą.

Ponieważ całka jest całką niewłaściwą to nie dla

wszystkich funkcji f(t) spełniających podane warunki jest ona zbieżna.

Całką Laplace’a funkcji f nazywamy całkę niewłaściwą postaci:

Transformacja Laplace’a c.d.

Funkcję f(t) dla których istnieje całka Laplace’a

nazywamy oryginałami natomiast odpowiadające i funkcje F(s) nazywamy transformatami Laplace’a.

Samą transformację Laplace’a zwaną także przekształceniem Laplace’a w środowisku inżynierskim często zapisujemy jako:

dt t

f e

t f L s

F ( ) { ( )}

( )









Transformacja odwrotna Laplace’a

Transformacja odwrotna Laplace’a dla klasy funkcji spełniających podane założenia wyraża się

wzorem:

ds e s i F

x f s

F

L ^st

iu x

iu

u x^



 



   lim ( )

2 ) 1 ( )}

(

1{



Transformacja Laplace’a – przykładowe funkcje

0 ),

( x x  f

a s

1

2

2 a

s a



2

2 a

s s



2

2 a

s a

 s

)}

( { )

(s L f x

F 

xn sn^!¹

n

eax

ax sin

ax cos

ax sinh

ax cosh

s

1

Transformata pochodnej:

Transformata całki:



 ^



0

) ( )}

(

{ f x e f x dx

L ^{n sx n}

) 0 ( ...

) 0 ( )

(s f ¹ s ¹ f

F

sⁿ  ^n   ^n



) 1 (

} )

( {

0

s s F

d f

L

t



• Linowość:

Własności transformacji Laplace’a

• Przesunięcie w dziedzinie zmiennej zespolonej:

)) (

( ))

( ( ))

( ( )}

( )

( )

(

{ af x bg x ch x aL f x bL g x cL h x

L     

gdzie a, b, c to współczynniki

) ( )}

(

{ f x F s

L  L { e

f ( x )}  F ( s  a )

• Zmiana skali

 

 



 

a F s

ax a f

L 1

)}

( {

Jeśli to

) ( )}

(

{ f x F s

L 

Jeśli to

Transformacja Laplace’a - przykład

Rozwiązywanie równania różniczkowego przy pomocy transformacji Laplace’a:

5 )

0 ( 2

3  y  e^ y 

dx

dy _x

 

1 ) 1

( 2 )]

0 ( )

( [ 3

2 3

 





 



 



  ^

s s Y y

s sY

e L dx y

L dy ^x

Krok 1: Transformacja obydwu stron równania różniczkowego

Uwzględniając warunek początkowy y(0) = 5 otrzymujemy:

) 1 ( 2 ] 5 ) ( [

3 sY s   Y s 

Transformacja Laplace’a – przykład c.d.

) 2 3

)(

1 (

16 ) 15

(

1 16 ) 15

( ) 2 3

(

1 15 ) 1

( ) 2 3

(





 



 



 





s s

s s Y

s s s

Y s

s s Y s

Krok 2: Wyrażanie Y(s) jako funkcji s

Zapisanie Y(s) w postaci ułamków prostych

2 3

18 1

) 1 (

; 18

; 2 1

3 1 )

2 3

)(

1 (

16 15

 



 





 

 

 





s s s

Y

B s A

B s

A s

s s

Transformacja Laplace’a – przykład c.d.

Krok 3: Odwrotna transformacja obydwu stron równania



 





 



 







 

 



 

 



 







666667 .

0 6 1

)} 1 ( {

666667 .

0 6 1

1 2

3 18 1

) 1 (

1 1

1

L s L s

s Y L

s s

s s s

Y

e at

a

L^ s   ^



 







1 1

x

x

e

e x

y ( )  

 6

Uwzględniając że:

Rozwiązanie równania wynosi:

Równanie Poissona

Równaniem Poissona nazywamy niejednorodne równanie Laplace’a

 

2 2 2

2 2

2

,

,

) , ( )

, (

z y

x

R z

y x

r

t r f

t r u



 



 



 











Równanie Poissona

Równanie Poissona możemy podać explicite dla przestrzeni 2 i 3 wymiarowej:

) , , ( )

, , ( )

, , ( )

, , (

) , ( )

, ( )

, (

2 2 2

2 2

2

2 2 2

2

z y x f z

y x z u

z y x y u

z y x x u

y x f y

x y u

y x x u

 

 



 





 

 





Równanie Poissona

Rozwiązanie równania Poissona wyrażamy za pomocą funkcji Greena:

' 4

) 1 ' , (

' )

' ( )

' , ( )

(

r r r

r G

dr r

f r

r G r

u

R

 

 



Rozwiązywanie równań różniczkowych zwyczajnych – metoda Eulera

Metoda Eulera pozwala na numeryczne rozwiązanie równania różniczkowego zwyczajnego rzędu

pierwszego postaci:

  ^{x , y , y}   ^{0 y}

dx f

dy  

Przykłady:

 

 

 

f

y y

dx e dy

y e

dx y dy

x x

x

2 3

. 1 ,

5 0

, 2 3

. 1

5 0

, 3

. 1 2























 

 

 

y y

y x y x

x f

e y

y x x

dx dy

y x

y dx x

e dy

2 2 2 2 2

2

) 3 sin(

, 2

5 0

) ,

3 sin(

2

5 0

), 3 sin(

2

 

 









Rozwiązywanie równań różniczkowych zwyczajnych – metoda Eulera

Dla x = 0 wartość y = y₀ przyjmując że x = x₀ = 0

y

Φ

krok h

x wartość prawdziwa y₁, wartość

przewidywana )

, (x₀ y₀

x1

Znając f(x, y) i mając dane wartości x₀ i y₀ z warunku początkowego y(x₀) = y₀ można obliczyć nachylenie funkcji f(x, y) do osi X w punkcie (x₀, y₀)

Rozwiązywanie równań różniczkowych zwyczajnych – metoda Eulera

) ,

(

tan _{0 0}

0 1

0

1 f x y

x x

y

y 



 









0

1

y f x , y x x

y   

 ^{x y}  ^h

f y

y





,

Po przekształceniu otrzymujemy:

Oznaczając x₁-x₀ jako krok h otrzymujemy:

Wykorzystując obliczaną wartość y₁ można obliczyć wartość y₂ dla x₂ jako:

 ^{x y}  ^h

f y

y





,

x

 x

 h

Rozwiązywanie równań różniczkowych zwyczajnych – metoda Eulera

 ^{x y}  ^h

f y

y





,

Można zatem wyprowadzić wzór rekurencyjny:

Φ

krok h

wartość prawdziwa

y_i+1, wartość przewidywana

y_i

x y

x_i x_i+1

Metoda Eulera - Przykład





^{ }

10 2067

.

2  ^{12 4}   ⁸ 



 ^  

 dt d

Kula nagrzana do temperatury 1200 K schładza się do temperatury 300K. Zakładając że w procesie

schładzania ciepło jest oddawane do otoczenia

jedynie przez radiację to temperatura kuli opisana jest równaniem:

Jaka jest temperatura kuli po 480 sekundach jej schładzania?

  ^{t ,    2 . 2067  10}

 ^

^{ 81  10}



f

Metoda Eulera – Przykład c.d.

Wzór rekurencyjny metody Eulera:

^

^{ }

^{ f}  ^t

^{, }

 ^h

Dla i = 0, t₀ = 0, Θ₀ = 1200:

Zakładamy krok h = 240

   

 

 

  ^K

f h

t f

09 . 106 240

5579 .

4 1200

240 10

81 1200

10 2067

. 2 1200

240 1200

, 0 1200

,

8 4

12 0

0 0

1











































240 240

0

0

1

    

 t t h

t

Metoda Eulera – Przykład c.d.

Dla i = 1, t₁ = 240, Θ₀ = 106.09:

   

 

 

 

f h

t f

32 . 110 240

017595 .

0 09

. 106

240 10

81 09

. 106 10

2067 .

2 09

. 106

240 09

. 106 ,

240 09

. 106 ,

8 4

12 1

1 1

2









































480 240

1 240

2     

 t t h t

Po wykonaniu dwóch iteracji metody Eulera otrzymujemy że temperatura kuli po 480 s wyniesie 110.32K

Metoda Eulera – Przykład c.d.

czas t(s)

dokładne rozwiązanie

temperatura Θ(K)

Porównanie rozwiązania dokładnego z rozwiązaniem otrzymanym przy użyciu metody Eulera.

Metoda Eulera – Przykład c.d.

Dobór odpowiedniego kroku h jest kluczowy dla odpowiedniej dokładności rozwiązania stosując metodę Eulera

 ⁴⁸⁰ 



rozmiar kroku h

480240 120 60 30

-987.81 110.32 546.77 614.97 632.77

252.54 82.964 15.566 5.0352 2.2864

czas t(s)

temperatura Θ(K)

dokładne rozwiązanie

Rozwiązywanie równań różniczkowych zwyczajnych – metoda Rungego - Kutty Metoda Rungego-Kutty pozwala na numeryczne rozwiązanie równania różniczkowego zwyczajnego pierwszego rzędu postaci:

 

 

dx f

dy  

Korzystając z rozwinięcie w szereg Taylora:

      ...

! 3

1

! 2

1 3

1 ,

3 2 3

1 ,

2 2 1

,

1   _   _   _  

 i i

y x i

i y x i

i y x i

i x x

dx y x d

dx x y x d

dx x y dy

y

i i i

i i i

   

 

! 3 ) 1

, (

! ' 2 ) 1

,

( ₁   ₁  ²  ₁  ³ 



 y_i f x_i y_i x_i x_i f x_i y_i x_i x_i f x_i y_i x_i x_i

Widać że metoda Eulera jest metodą Rngego-Kutty pierwszego rzędu

 ^{x y}  ^h

f y

y





,

Rozwiązywanie równań różniczkowych zwyczajnych – metoda Rungego - Kutty Wzór dla metody Rungego-Kutty drugiego rzędu będzie wyglądał następująco:

   

1 ,

! 2

, y h 1 f x y h x

f y

y_i  _i  _{i i}   _{i i}

 

 

 

 

1 ,

! 4 , 1

! 3 , 1

! 2

, y h 1 f x y h f x y h f x y h x

f y

y_i  _i  _{i i}  _{i i}  _{i i}  _{i i}

dla metody czwartego rzędu:

Jak wyznaczyć pochodne f’ metody stopnia drugiego i f’, f’’, f’’’ dla metody stopnia czwartego?

Rozwiązywanie równań różniczkowych zwyczajnych – metoda Rungego - Kutty

 ^{a k a k}  ^h

y

y







   

1 ,

! 2

, y h 1 f x y h x

f y

y_i  _i  _{i i}   _{i i}

Dla metody Rungego-Kutty drugiego rzędu wzór:

można zapisać jako:





f

k₁  ,

gdzie: ^k₂  ^f





aby wyznaczyć współczynniki a₁, a₂, p₁, q₁₁ należy rozwiązać kład równań:

2

1

1

 a 

a 2

1

1 2

p 

a 2

1

11 2

q  a

zazwyczaj wartość a₂ wybiera się aby rozwiązać pozostałe

Rozwiązywanie równań różniczkowych zwyczajnych – metoda Rungego - Kutty Zazwyczaj a₂ przyjmuje jedną z trzech wartości: ½, 1, 2/3

Metoda Heun’a Metoda midpoint Metoda Raltson’a 2

1

2  a

2 1

1 

a p₁ 1 q₁₁ 1

h k k

y

y_{i i} 



 



 



1  1 2

2 1 2

1





f

k₁  ,





f

k₂  _i  , _i  ₁

2 1 a

1  0

a 2

1

1 

p 2

1

11  q

h k y

y_i1  _i  ₂





f

k₁  ,



 



  

 f x h y k h

k_{2 i i 1}

2 , 1

2 1

3 2

2  a

3 1

1 

a 4

3

1 

p 4

3

11  q

h k k

y

y_{i i} 



 



 



1  1 2

3 2 3

1





f

k₁  ,



 



  

 f x h y k h

k_{2 i i 1}

4 , 3

4 3

Metoda Rungego - Kutty - Przykład





^{ }

10 2067

.

2  ^{12 4}   ⁸ 



 ^  

 dt d

Kula nagrzana do temperatury 1200 K i schładza się do temperatury 300K. Zakładając że w procesie

schładzania ciepło jest oddawane do otoczenia

jedynie przez radiację to temperatura kuli opisana jest równaniem:

Jaka jest temperatura kuli po 480 sekundach jej schładzania?

  ^{t ,    2 . 2067  10}

 ^

^{ 81  10}



f

Metoda Rungego - Kutty – Przykład c.d.

 





f k

t f k

h k k

i i

i i i i

1 2

1

2 1

1

, ,

2 1 2

1











 



 



 







Dla metody Heun’a: 



 







1  f t  f    ^    

k _o

       





10 2067

. 2

09 . 106 , 240 240

5579 .

4 1200

, 240 0

,

8 4

12 1 0

0 2































f f

h k h

t f

k 

Dla i = 0, t₀ = 0, Θ₀ = Θ(0) = 1200:

   

 

h k

k 0.017595 240

2 5579 1

. 2 4

1200 1 2

1 2

1

2 1

0

1 



 



  



 



 



 





