Numerikus módszerek 5. Közönséges differenciálegyenletek numerikus megoldása

Numerikus módszerek 5. Közönséges differenciálegyenletek numerikus megoldása

Kezdeti és peremérték feladatok Az Euler-módszer Az Euler-módszer javításai Runge-Kutta-módszerek Lineáris többlépéses módszerek Peremérték feladatok másodrendű differenciálegyenletekre

1

Kezdeti érték feladatok

y( x)  f ( x, y( x)), y( x0 )  y0

( x0  x  x0  A) (kezdeti feltétel)

(ahol f adott, kétváltozós függvény, mely kielégíti a Lipschitz-feltételt). Példa (oldat hígítása). Legyen egy V térfogatú tartály valamilyen vízben oldott anyaggal (pl. sóval). Kezdjünk el folyamatosan, egyenletesen tiszta vizet tölteni a tartályba felül, és hagyjuk a sós vizet alul ugyanilyen ütemben távozni. A tartályban nyilván egyre hígabb lesz a sóoldat. A hígulás folyamatát akkor az alábbi differenciálegyenlet írja le:

c(t )  

Q c(t ) V

ahol c(t ) a t időpillanatban érvényes sókoncentráció a tartályban, és Q a hozzáfolyás hozama, (időegység alatt a tartályba beáramló víz mennyisége (térfogata)). Kezdeti feltétel: c(0)  c0 ahol c0 a kezdeti (a t  0 időpillanatban érvényes) sókoncentráció.

2

Peremérték feladatok

y( x)  f ( x, y( x), y( x)), y( x0 )  y0 , y( x1 )  y1

( x0  x  x1 ) (peremfeltételek)

(ahol f adott, háromváltozós függvény, mely kielégíti a Lipschitz-feltételt). Peremfeltételként a deriváltértékek is előírhatók a perempontokban, vagy az érték és a derivált adott lineáris kombinációja is. Példa (elektromos áram vezetőben). Tekintsünk egy vékony, hosszú, de nem feltétlen mindenütt ugyanolyan vastagságú és/vagy anyageloszlású elektromos vezetékdarabot. Kapcsoljunk áramforrást a vezeték két végpontjára. Az elektromos feszültség eloszlását a vezeték hossza mentén ekkor az alábbi differenciálegyenlet írja le:

(  (U )  0 ahol U (x) az elektromos potenciál a vezeték x pontjában,  (x) a relatív vezetőképesség ugyanitt. A peremfeltételek pl: U ( x0 )  U 0 , U ( x1 )  U1 ahol U 0 , U1 a potenciál a vezeték két végpontjában. 3

Kezdeti érték feladatok megoldási elvei Hasznos segédeszközök: Taylor-formula egyváltozós függvényekre:  f ( k ) ( x) f ( x) f ( x) 2 f ( x) 3 f ( x  h)  f ( x )  h h  h  ...   hk 1! 2! 3! k! k 0 Taylor-formula kétváltozós függvényekre:

1  f ( x, y ) f ( x , y )  1   2 f ( x , y ) 2  2 f ( x, y )  2 f ( x, y ) 2  f ( x  h, y   )  f ( x , y )   h   h 2 h    2 2  x  y  x  y   y 1! 2! x  3 1   3 f ( x, y ) 3  3 f ( x, y ) 2  3 f ( x, y ) 2  f ( x, y ) 3    h 3 2 h  3 h     ... 3 2 3 x y xy y 3! x 

1   k   k f ( x, y ) k  k   k f ( x, y ) k 1 k   k f ( x, y ) k      h    k 1 h   ...       1  x y  k  y k k!  0  x k  

1 k  k   k f ( x, y ) k  j      k j j h  j y k  0 k! j  0   x

j

4

Számítási rács: xk : x0  kh (k  0,1,..., N ) , ahol h a lépésköz. Cél: az y megoldás értékeinek közelítése a rácspontokban. 1. A deriváltak közelítése véges differenciákkal Elsőrendű derivált: y y  Előrelépő séma: y( xk ) ~ k 1 k Hibája: O(h) , h mert a Taylor-formulából: yk 1  y( xk 1 )  y( xk  h)  yk  y( xk )h  O(h 2 ) y  yk 1  Visszalépő séma: y( xk ) ~ k Hibája: O(h) , h mert a Taylor-formulából: yk 1  y( xk 1 )  y( xk  h)  yk  y( xk )h  O(h 2 ) y  yk 1  Centrális séma: y( xk ) ~ k 1 Hibája: O(h 2 ) , 2h 1 mert a Taylor-formulából yk 1  y ( xk 1 )  y ( xk  h)  yk  y( xk )h  y( xk )h 2  O(h3 ) 2 1 és hasonlóan: yk 1  y ( xk 1 )  y ( xk  h)  yk  y( xk )h  y( xk )h 2  O(h3 ) 2

5

Másodrendű derivált:  Centrális séma: y( xk ) ~

yk 1  2 yk  yk 1 h2

Hibája: O(h 2 ) ,

mert a Taylor-formulából:

1 1 y( xk )h 2  y( xk )h3  O(h 4 ) 2 6 1 1 yk 1  y ( xk 1 )  y ( xk  h)  yk  y( xk )h  y( xk )h 2  y( xk )h3  O(h 4 ) 2 6 Összeadva a két egyenlőséget, az állítás adódik. yk 1  y ( xk 1 )  y ( xk  h)  yk  y( xk )h 

2. A differenciálegyenlet integrálásával y ( xk 1 )  y ( xk ) 

x k 1

 f ( x, y( x)) dx

xk

és a jobb oldali integrált egy közelítő integrálformulával helyettesítjük. h Példa (trapézformula alkalmazása): yk 1  yk   ( f ( xk , yk )  f ( xk 1, yk 1)) 2

6



7

Az Euler-módszer Modellfeladat: y( x)  f ( x, y( x)), y( x0 )  y0

( x0  x  x0  A) (kezdeti feltétel)

Számítási rács: ekvidisztáns, h lépésközű: x0 , x1, x2 ,..., x N . A megoldás deriváltját a rácspontokban elsőrendű véges differenciákkal közelítjük. Explicit séma (előrelépő séma használata): y  yk y ( xk )  k 1 : f ( xk , yk )  yk 1 : yk  h  f ( xk , yk ) h

(k  0,1,...,N  1)

Implicit séma (visszalépő séma használata): y  yk y ( xk 1 )  k 1 : f ( xk 1 , yk 1 )  yk 1 : yk  h  f ( xk 1, yk 1 ) h

(k  0,1,...,N  1)

8

Konzisztencia, stabilitás, konvergencia

y ( xi 1 )  y ( xi )  f ( xi , y ( xi )) (i  0,1,..., N  1) h Egy módszer p-edrendben konzisztens, ha gi  O(h p ) (i-től függetlenül). A lokális hibatagok: gi :

A globális hibatagok: ei : y( xi )  yi

(i  0,1,..., N  1)

Egy módszer p-edrendben konvergens, ha ei  O(h p ) (i-től függetlenül). Azt mondjuk, hogy egy módszer stabil, ha a globális hiba a lokális hibával becsülhető felülről: i 1    | ei |  C  | e0 |  h  | g j |  (i  0,1,...,N  1)   j 0   Ha egy módszer p-edrendben konzisztens, és stabil, akkor p-edrendben konvergens is. Az Euler-módszer elsőrendben konvergens.

9



10

Az Euler-módszer javításai Az eredeti differenciálegyenletet átírjuk a következő alakba: y ( xk 1 )  y ( xk ) 

x k 1

 f ( x, y( x)) dx

xk

b

Trapézformula: Használva az  F ( x)dx  (b  a)  a

F (a)  F (b) integrálközelítő formulát: 2

h yk 1 : yk    f ( xk , yk )  f ( xk 1, yk  hf ( xk , yk ))  2 Másképp felírva: yk* 1 : yk  h  f ( xk , yk )



h yk 1 : yk   f ( xk , yk )  f ( xk 1, yk* 1 2



A módszer másodrendben konzisztens.

11

a b Érintőformula: Használva az  F ( x)dx  (b  a)  F   integrálközelítő formulát: 2   a b

yk 1 : yk  h  f ( xk 1/ 2 , yk 

h f ( xk , yk )) 2

Másképp felírva:

h yk 1/ 2 : yk   f ( xk , yk ) 2 yk 1 : yk  h  f ( xk 1/ 2 , yk 1/ 2 ) A módszer másodrendben konzisztens.

12

Az aszimptotikus stabilitás öröklődése Tekintsük az

y   Ay

modellfeladatot ( A  0 ), melynek azonosan 0 megoldása aszimptotikusan stabil (tetszőleges y megoldásra teljesül, hogy lim y ( x)  0 , mivel y( x)  y(0)  e  Ax ) x  

Legyen most h rögzített, és vizsgáljuk, hogy a közelítő megoldásra yk  0 mikor teljesül. Explicit Euler-módszer:

yk 1 : yk  Ahyk

Az aszimptotikus stabilitás nem minden h mellett öröklődik (feltételes stabilitás), mert yk  yk 1  Ahyk 1  (1  Ah) yk 1  (1  Ah) 2 yk  2  ...  (1  Ah) k y0 .

Így yk  0 akkor teljesül minden y0 mellett, ha | 1  Ah | 1, azaz  1  1  Ah  1. Az aszimptotikus stabilitás tehát akkor öröklődik, ha

0h

2 A 13

Implicit Euler-módszer:

yk 1 : yk  Ahyk 1

Az aszimptotikus stabilitás minden pozitív h mellett öröklődik (feltétel nélküli stabilitás), mert yk : yk 1  Ahyk miatt:

yk 

1 1 1 yk 1  y  ...  y0 k 2 2 k 1  Ah (1  Ah) (1  Ah)

Így yk  0 minden (pozitív) h mellett teljesül.

14

A pozitivitás megőrzése Tekintsük az

y   Ay

modellfeladatot ( A  0 ), melynek minden pozitív kezdeti feltételhez tartozó megoldása mindenütt pozitív (mert y( x)  y(0)  e  Ax ). Legyen most h rögzített, és vizsgáljuk, hogy a közelítő megoldásra yk  0 mikor teljesül, ha y0  0 . Explicit Euler-módszer:

yk 1 : yk  Ahyk

A pozitivitás nem őrződik meg minden h mellett, mert yk  yk 1  Ahyk 1  (1  Ah) yk 1  (1  Ah) 2 yk  2  ...  (1  Ah) k y0 .

Így yk  0 akkor teljesül minden y0  0 mellett, ha 0  1  Ah  1. A pozitivitás tehát akkor őrződik meg, ha 1 0h A 15

Implicit Euler-módszer:

yk 1 : yk  Ahyk 1

A pozitivitás minden pozitív h mellett megőrződik, mert yk : yk 1  Ahyk miatt:

yk 

1 1 1 yk 1  y  ...  y 2 k 2 k 0 1  Ah (1  Ah) (1  Ah)

Így yk  0 minden (pozitív) h mellett teljesül.

16



17

Runge-Kutta-módszerek k1 : f ( xi , yi ) k 2 : f ( xi  h  a2 , yi  h  b21k1 ) k3 : f ( xi  h  a3 , yi  h  b31k1  h  b32 k 2 ) ... k s : f ( xi  h  as , yi  h  bs1k1  h  bs 2 k 2  ...  h  bs , s 1k s 1 ) yi 1 : yi  h  (c1k1  c2 k 2  ...  cs k s )

Egy-egy konkrét módszert jellemző paraméterek: s a2 , a3 ,...,as b21; b31, b32 ; b41, b42 , b43 ; ... bs1, bs 2 ,...,bs , s 1 c1, c2 ,...,cs

y ( xi 1 )  y ( xi ) s   c jk j A paraméterek úgy választandók meg, hogy a képlethiba tagjai, a h j 1 számok O(h p ) nagyságrendűek legyenek, ahol p adott szám (a módszer rendje). 18

Speciális Runge-Kutta-módszerek Elsőrendű Runge-Kutta-módszer: megegyezik az Euler-módszerrel. Másodrendű Runge-Kutta-módszerek: az érintő- ill. a trapézformulával javított Euler-módszerek. Egy harmadrendű Runge-Kutta-módszer: h h k1 : f ( xi , yi ), k 2 : f ( xi  , yi  k1 ), k3 : f ( xi  h, yi  h  k1  2h  k 2 ) 2 2 h yi 1 : yi   (k1  4k 2  k3 ) 6 Egy negyedrendű Runge-Kutta-módszer:

k1 : f ( xi , yi ),

h h k 2 : f ( xi  , yi  k1 ), 2 2

h h k3 : f ( xi  , yi  k 2 ), k 4 : f ( xi  h, yi  h  k3 ) 2 2 h yi 1 : yi   (k1  2k 2  2k3  k 4 ) 6 19

Runge-Kutta-módszerek rendje Adott s paraméter mellett az elérhető legmagasabb rend: ( s  1,2,3,4)  s  s  1 ( s  5,6,7)  p ( s  8,9)  s  2   s  3 ( s  10)

20



21

Lineáris többlépéses módszerek A lineáris k-lépéses módszerek általános alakja: k

k

j 0

j 0

  j yi  j  h   j f ( xi  j , yi  j )

(i  0,1,..., N  k )

y0 a kezdeti feltételből adott; y1, y2 ,..., yk 1 meghatározása teljesen különálló módszerrel kell, hogy történjék (pl. egy Runge-Kutta-módszer segítségével). Ez az induló eljárás feladata. Egy-egy konkrét módszert jellemző paraméterek:

k ,  0 , 1,..., k ,  0 , 1,...,  k

(ahol  k  0 ),

melyekkel definiáljuk a módszerre jellemző polinomokat: k

 ( z ) :   j z , j 0

j

k

 ( z ) :   j z j j 0

22

Konzisztencia A paraméterek úgy választandók meg, hogy a képlethiba tagjai, azaz a k 1 k g i :   j y ( xi  j )    j f ( xi  j , y ( xi  j )) h j 0 j 0 számok O(h p ) nagyságrendűek legyenek, ahol p adott szám (a módszer rendje). A módszer rendje (legalább) p, ha a következő feltételek valamelyike teljesül: (a) C0  C1  C2  ...  C p  0 , k

1 k 1 k 1 k 2 1 k ahol C0 :   j , C1 :  j j    j , C2 :  j  j   j j ,…, 1! j  0 0! j  0 2! j  0 1! j  0 j 0 k 1 k p 1 p 1 C p :  j j   j  j (hibakonstansok) p! j  0 ( p  1)! j  0 (b) a  ( z ) /  ( z ) racionális törtfüggvény a z = 1 hely körül a (komplex változós) logaritmus ( z)  log z  O(( z  1) p 1 ) . függvényt (p +1)-edrendben közelíti, azaz  ( z)

23

Stabilitás és konvergencia A módszer stabil, ha a  polinom kielégíti a gyökkritériumot, azaz ha minden gyökének abszolút értéke legfeljebb 1, és minden 1 abszolút értékű gyöke csak egyszeres gyök. Ekkor, amennyiben a módszer rendje p, a módszer p-edrendben konvergens is. Ha a  polinom kielégíti a gyökkritériumot, akkor az elérhető legmagasabb rend: p  k  2 ha k páros, ill. p  k  1, ha k páratlan. (Dahlquist tétele)

24

Adams-típusú lineáris k-lépéses módszerek

 ( z )  z k  z k 1 Explicit, k-adrendű módszerek (Adams-Bashforth-módszerek):

k  1:

yi 1 : yi  h  f ( xi , yi )

k  2:

h yi 1 : yi  (3 f ( xi , yi )  f ( xi 1, yi 1 )) 2

k  3:

yi 1 : yi 

k  4: yi 1 : yi 

h (23 f ( xi , yi )  16 f ( xi 1, yi 1 )  5 f ( xi  2 , yi  2 )) 12

h (55 f ( xi , yi )  59 f ( xi 1, yi 1 )  37 f ( xi  2 , yi  2 )  9 f ( xi 3 , yi 3 )) 24

25

Adams-típusú lineáris k-lépéses módszerek

 ( z )  z k  z k 1 Implicit, (k  1) -edrendű módszerek (Adams-Moulton-módszerek):

k  1: h yi 1 : yi  ( f ( xi 1, yi 1 )  f ( xi , yi )) 2

k  2:

yi 1 : yi 

k  3: yi 1 : yi 

h (5 f ( xi 1, yi 1 )  8 f ( xi , yi )  f ( xi 1, yi 1 )) 12

h (9 f ( xi 1, yi 1 )  19 f ( xi , yi )  5 f ( xi 1, yi 1 )  f ( xi  2 , yi  2 )) 24

26



27

Peremérték feladatok másodrendű differenciálegyenletekre Modellfeladat: Legyen  mindenütt pozitív, q mindenütt nemnegatív adott függvény:  ( ( x)u ( x))  q( x)u ( x)  f ( x), u (0)  a, u (1)  b

(0  x  1)

Számítási rács: xk : kh (k  0,1,..., N ) , ahol h a lépésköz: h : 1 / N . A másodrendű tag közelítése centrális sémával történik: 1 u u u  uk 1  1 (u)( xk ) ~  k 1 k 1 k  k k   2 k 1uk 1  (k 1  k )uk  k uk 1  h h h  h

u0 , u N adottak; a többi u k egy lineáris egyenletrendszer megoldásával kapható:  k uk 1  (k 1  k  h2qk )uk  k 1uk 1  h2 f k

(k  1,2,..., N  1)

Az egyenletrendszer mátrixa 3-átlós (sőt, ha mindegyik qk pozitív, akkor diagonálisan domináns is), így a megoldás gazdaságosan elvégezhető (O(N) műveletigénnyel).

28

Numerikus módszerek 5. Közönséges differenciálegyenletek numerikus megoldása

Recommend Documents