Matematika pro chemické inženýry

(1)

Drahoslava Janovsk ´a

7. Numerick ´e ˇre ˇsen´ı ODR - po ˇc ´ate ˇcn´ı ´ uloha:

Eulerova metoda, Rungova-Kuttova metoda 2. a 4. ˇr ´adu.

(2)

(ˇz ´adn ´e oznaˇcen´ı)

F Pˇr´ıklady k procviˇcen´ı - dobrovoln ´e

F Pro studenty, kteˇr´ı cht ˇej´ı v ˇed ˇet v´ıc. Tato l átka se nebude pˇredn áˇset, nebude v p´ısemk ách, nebude se zkouˇset.

(3)

Obsah

1 Jednokrokov ´e metody Eulerova metoda

Rungeovy–Kuttovy metody

2 V´ıcekrokov ´e metody

Line ´arn´ı v´ıcekrokov ´e metody

3 Stabilita k −krokov ´ych metod

4 Literatura k dal ˇs´ımu studiu

(4)

Numerick é ˇreˇsen´ı diferenci áln´ıch rovnic – nezbytnost pro inˇzenýrsk é aplikace Budeme uvaˇzovat jednu dif. rovnici i ˇreˇsen´ı soustav dif. rovnic. Dif. rovnice vyˇs´ıch ˇr ád ˚u lze pˇrev ést na syst ém rovnic 1. ˇr ádu. Budeme se zabývat zejm éna rovnicemi 1. ˇr ádu, a to úlohou poˇc áteˇcn´ı i okrajovou.

(5)

Jednokrokov ´e metody pro ˇre ˇsen´ı po ˇc ´ate ˇcn´ı ´ ulohy

Reˇs´ıme poˇc ´ateˇcn´ı (Cauchyovu) ´ulohuˇ

y⁰=f (x , y ), y (x0) =y0.

Jednoznaˇcnost a existence ˇreˇsen´ı:

Je-li f (x , y ) spojit ´a v Ω = {(x , y ), |x − x0| ≤ a, a > 0, |y − y0| ≤ b, b > 0} a oznaˇc´ıme-li

M = max

(x ,y )∈Ω|f (x, y )| , α =min(a, b M) , pakexistuje ˇreˇsen´ırovnice y⁰=f (x , y ) definovan ´e v intervalu (x0− α, x0+ α).

Je-li nav´ıc f Lipschitzovsk ´a, t.j.

|f (x, y1) −f (x , y2)| ≤L |y1− y2| L > 0, ∀(x, y1), (x , y2) ∈ Ω , jeˇreˇsen´ı urˇceno jednoznaˇcn ˇe.

Pozn ´amka: L . . . Lipschitzovsk ´a konstanta.

(6)

Taylor ˚ uv rozvoj ˇre ˇsen´ı

Reˇs´ıme ´ulohu, tj. hled ´ame ˇreˇsen´ı y = y (t) tak, abyˇ dy

dt =y⁰=f (t, y ), y (0) = y0, 0 < t < T .

Necht’ uˇz jsme spoˇc´ıtali ˇreˇsen´ı yn:=y (tn)v ˇcase tna chceme naj´ıt ˇreˇsen´ı v ˇcase tn+1. Necht’ h := tn+1− tnje pˇr´ısluˇsn´y ˇcasov´y krok.

Taylor ˚uv rozvoj ˇreˇsen´ı:yn+1v bod ˇe yn: yn+1=yn+h · yn⁰

|{z}

+h² 2yn⁰⁰+h³

6yn⁰⁰⁰+ . . .

=f (tn,yn)z dif. rovnice Co yn⁰⁰?

y⁰⁰= dy dt

0

=∂f

∂t · 1 + ∂f

∂y · y⁰=ft+fy· f y⁰⁰⁰= ∂

∂t(ft+f fy) + ∂

∂y(ft+f fy)f = ftt+2f fyt+ftfy+f fy²+f²fyy, . . . atd.

(7)

Diskretiza ˇcn´ı chyba, ˇr ´ad metody

Jak se pˇribliˇzn é ˇreˇsen´ı z´ıskan é danou metodou liˇs´ı od pˇresn ého ˇreˇsen´ı?Tedy jak á jeglob áln´ı diskretizaˇcn´ı chyba

ei=y (xi) −yi ?

Pˇri numerick ém ˇreˇsen´ı diferenci áln´ı rovnice se v kaˇzd ém kroku dopouˇst´ıme lok áln´ı diskretizaˇcn´ı chyby. Glob áln´ı diskretizaˇcn´ı chyba je výsledkem nakupen´ı lok áln´ıch chyb. O metod ách, u kterých se daˇr´ı udrˇzet chybu malou vzhledem k pˇresn ému ˇreˇsen´ı, mluv´ıme jako ostabiln´ıch metod ách.

Pro popis rychlosti konvergence pouˇz´ıv áme pojemˇr ád metody. ˇR ád metody je pˇrirozen é ˇc´ıslo p, takov é, ˇze pro mal á h je lok áln´ı diskretizaˇcn´ı chyba ˇr ádov ˇe velikosti h^p+1.

Napˇr´ıklad

yn+1=yn+hy_n⁰+ O(h²)

| {z }

znamen ´a, ˇze lim

h→O₊

yn+1

h² =k 6= 0, lok ´aln´ı diskretizaˇcn´ı chyba

metoda je ˇr ádu p = 1. N ásleduj´ıc´ıEulerova metodaje metoda prvn´ıho ˇr ádu.

(8)

Eulerova metoda

Eulerova metoda

Na prvn´ıch dvou ˇclenech Taylorova rozvoje je zaloˇzenaEulerova metoda:

Zvol´ımy (0) := y0a konstruuji posloupnost yn+1=yn+hf (tn,yn),

t.j. v bod ˇe (tn,yn)nahrad´ım funkci y teˇcnou k integr ´aln´ı kˇrivce v bod ˇe (tn,yn).

Pˇresnost v´ypoˇctu:Eulerova metoda je metoda 1. ˇr ´adu.

Pozn ´amka Krok metody h lze m ˇenit s jednotliv´ymi iteracemi (adaptivn´ı volba kroku):

yn+1=yn+hnf (tn,yn) .

(9)

Eulerova metoda

F Pˇr´ıklad

Pˇr´ıklad Najd ˇete Eulerovou metodou ˇreˇsen´ı n ásleduj´ıc´ı poˇc áteˇcn´ı úlohy v bod ˇe (ˇcase) t = 3,

y⁰=0, 3 y sin(t) , y (1) = 2 . Pro jednoduchost uvaˇzujme n = 4, je tedy h = 0.5 .

yj+1=yj+h f (tj,yj), kde f (tj,yj) =0, 3yjsin(tj) .

tj krok 0.5 krok 0.005

1 2 2

1.5 2.252441295 2.30249902026881692 2 2.589461130 2.66460601831410714 2.5 2.942649681 2.99089235783755570 3 3.206813761 3.16533517440834976

(10)

Eulerova metoda

F Pˇr´ıklad

Eulerovou metodou ˇreˇste diferenci ´aln´ı rovnici y⁰=t − 2y , y (0) = 1 . Pˇresn ´e ˇreˇsen´ı: y (t) =1

4 h

2t − 1 + 5e^−2t i

.Integraˇcn´ı krok h = 0.2, tˇri iterace.

Eulerova metoda pˇresn ´e ˇr. chyba j t_j f (t_j−1,y_j−1) y_j=y_j−1+hf (t_j−1,y_j−1) y (t_j) y_j− y (tj)

0 0.0 poˇc.podm. = 1.0000 0 0

1 0.2 0 − 2 · 1 = −2.000 1.0 + (0.2)(−2.0) = 0.6000 0.6879 −0.0879 2 0.4 0.2 − (2)(0.6) = −1.000 0.6 + (0.2)(−1.0) = 0.4000 0.5117 −0.1117 3 0.6 0.4 − (2)(0.4) = −0.400 0.4 + (0.2)(−0.4) = 0.3200 0.4265 −0.1065

(11)

Eulerova metoda

F

Na obr ázku je srovn án´ı pˇresn ého ˇreˇsen´ı a numerick ého ˇreˇsen´ı (•) Eulerovou metodou. Integr áln´ı kˇrivky z(t) startuj´ı vˇzdy z bod ˚u numerick ého ˇreˇsen´ı jako z nov é poˇc áteˇcn´ı podm´ınky pro danou rovnici.

(12)

Eulerova metoda

F Stabilita Eulerovy metody

Reˇsme modelovou ´ulohuˇ

y⁰= λy , λkonstanta (1)

Eulerova metoda =⇒

yn+1=yn+ λh yn, t.j. yn+1=yn(1 + λ h) . Tedy

yn=yn−1(1 + λ h) = yn−2(1 + λ h)²= . . . =y0(1 + λ h)ⁿ. Pro λ = λ1+ iλ2 . . . komplexn´ı, dostaneme

yn=y0(1 + λ1h + iλ2h

| {z }

)ⁿ=y0σⁿ.

σ . . . tzv. amplifikaˇcn´ı faktor

Numerick é ˇreˇsen´ı jestabiln´ı(t.j. z ˚ustane omezen é i pro rostouc´ı n (velk é)), jestliˇze|σ| ≤ 1.

(13)

Eulerova metoda

F

Necht’ v rovnici (1) je

λ = λ1+ iλ2, λ1≤ 0 , σ =1 + λ1h + iλ2h .

Pakoblast stability pro Eulerovu metoduje ˇc ´ast lev ´e poloroviny komplexn´ı roviny, a tovnitˇrek kruhu

|σ|²= (1 + λ1h)²+ λ²2h²=1 .

Pro libovolnou hodnotu λh v lev é polorovin ˇe komplexn´ı roviny vn ˇe tohoto kruhu je numerick é ˇreˇsen´ı rostouc´ı nade vˇsechny meze (blowing up), zat´ımco pˇresn é ˇreˇsen´ı kles á. Chceme-li m´ıt ˇreˇsen´ı stabiln´ı, mus´ımezmenˇsit htak, aby λh leˇzelo v tomto kruhu.

(14)

Eulerova metoda

F Implicitn´ı (zp ˇetn ´a) Eulerova metoda

V rovnici se vyskytuje yn+1implicitn ˇe:

yn+1=yn+h f (tn+1,yn+1)

Nevýhoda: metoda je výpoˇcetn ˇe n ároˇcn ˇejˇs´ı neˇz explicitn´ı Eulerova metoda Výhoda: je stabiln ˇejˇs´ı, n ˇekdy lze s výhodou vyuˇz´ıt linearizaci f

Pˇr´ıklad Aplikujme na modelovou ´ulohu (1) implicitn´ı Eulerovu metodu:

yn+1=yn+ λhyn+1 =⇒ yn+1= 1

1 − λhyn. Tedy yn= 1

1 − λhyn−1=

1

1 − λh

2

yn−2= . . . =

1

1 − λh

n

y0. Dostaneme:

yn= σⁿy0, σ = 1 1 − λh.

(15)

Eulerova metoda

F θ–metody

Lze definovat n ´asleduj´ıc´ı jednoparametrickou tˇr´ıdu jednokrokov´ych metod, tzv.θ−metody:

Pro danou poˇc ´ateˇcn´ı aproximaci y0definujeme yn+1jako konvexn´ı kombinaci f (tn,yn)a f (tn+1,yn+1)(θ . . .parametr):

yn+1=yn+h[(1 − θ)f (tn,yn) + θf (tn+1,yn+1)] ,n = 0, 1, . . . , N − 1 , θ ∈ h0, 1i .

• θ =0 =⇒ yn+1=yn+hf (tn,yn) . . . (explicitn´ı) Eulerova metoda

• θ =1 =⇒ yn+1=yn+hf (tn+1,yn+1) . . . implicitn´ı Eulerova metoda

• θ = 1

2 =⇒ yn+1=yn+1

2h[f (tn,yn) +f (tn+1,yn+1)]

. . . lichob ˇeˇzn´ıkov ´e pravidlo.

Lze uk ´azat, ˇze θ−metoda je explicitn´ı pro θ = 0 a je implicitn´ı pro 0 < θ ≤ 1 .

(16)

Rungeovy–Kuttovy metody

Rungeovy–Kuttovy metody

Rungeovy–Kuttovy metody (RK)– pˇresn ˇejˇs´ı neˇz Eulerovy metody:

explicitn´ı: ˇreˇsen´ı v ˇcase tn+1vypoˇcteme z hodnot yn,f (tn,yn)a z f (t, y ) vyˇc´ıslen ´e v bod ˇe leˇz´ıc´ım mezi body tna tn+1

=⇒ lepˇs´ı pˇresnost, protoˇze vyuˇzijeme v´ıce informac´ı o funkci f . implicitn´ı: obvykle vedou na ˇreˇsen´ı neline árn´ıch algebraických rovnic, ale pracnost je vyv áˇzena lepˇs´ı numerickou stabilitou.

(17)

RK metody 2. ˇr ´adu

Reˇsme rovnici yˇ ⁰=f (t, y ) .

V ˇcasov ´em kroku tn+1dostaneme ˇreˇsen´ı z rovnice

yn+1=yn+ γ1k1+ γ2k2, (2) kde

k1 = hf (tn,yn)

k2 = hf (tn+ αh, yn+ βk1) , α, β, γ1, γ2∈ R .

Konstanty α, β, γ1, γ2 mus´ıme urˇcit tak, aby metoda m ˇela co nejvyˇsˇs´ı ˇr ád pˇresnosti (byla co nejvyˇsˇs´ıho ˇr ádu). Abychom zjistili ˇr ád pˇresnosti, vyuˇzijeme Taylor ˚uv rozvoj y (tn+1):

yn+1=yn+h yn⁰

|{z}

+h² 2 yn⁰⁰

|{z}

+ · · · =⇒

f (tn,yn) ft+f fy

yn+1=yn+hf (tn,yn) +h²

2(ftn+fnfyn)+· · · (3) a porovn ´ame koeficienty v (2) a (3) .

(18)

Taylorova ˇrada pro funkci dvou prom ˇenn´ych k2=hf (tn+ αh, yn+ βk1) =⇒

k2=h

f (tn,yn) + αhft_n+ βk1fy_n+ O(h²) . Pozn ´amka Symbol O (velk ´e O)

g(h) = O(h^p) ⇐⇒ |g(h)| ≤ C · h^p, C konstanta nez ávisl á na h . Dostaneme yn+1=yn+ (γ1+ γ2)hfn+ γ2βh²fnfy_n+ γ2αh²ft_n+ · · · Porovn áme s (3) a dostaneme tˇri neline árn´ı rovnice pro 4 nezn ám é:

γ1+ γ2=1 , γ2α = 1

2, γ2β = 1 2. Necht’ α ∈ R je parametr. Pak γ2= 1

2α, β = α , γ1= (1 − 1 2α) .

(19)

Dostaneme Rungeovy-Kuttovy metody 2. ˇr ´adu:

k1 = hf (tn,yn)

k2 = hf (tn+ αh, yn+ αk1) yn+1 = yn+ (1 − 1

2α)k1+ 1 2αk2. Zvol´ım α a dostanu metodu. Napˇr.

α = 1

2 =⇒ γ2=1 , β = 1

2, γ1=0 =⇒

yn+1=yn+k2=yn+hf (tn+1

2h, yn+1 2k1) .

Pozn ánka RK metoda 2. ˇr ádu vyˇzaduje v kaˇzd ém kroku dvakr át vyˇc´ıslen´ı funkˇcn´ıch hodnot.

(20)

RK metody 4. ˇr ´adu

Pro ˇreˇsen´ı poˇc áteˇcn´ı úlohy se nejv´ıce pouˇz´ıv áRK metoda 4. ˇr ádu:

yn+1=yn+1 6k1+1

3(k2+k3) +1 6k4, kde

k1 = hf (tn,yn)

k2 = hf (tn+^h₂,yn+¹₂k1) k3 = hf (tn+^h₂,yn+¹₂k2) k4 = hf (tn+h, yn+k3) . V kaˇzd ´em kroku potˇrebujeme 4kr ´at vyˇc´ıslit funkˇcn´ı hodnoty.

Aˇckoliv pracn á, je RK metoda 4. ˇr ádu stabiln´ı a velmi pˇresn á. Je tak é snadno programovateln á, protoˇze nevyˇzaduje ˇz ádn é derivov án´ı, ale jen výpoˇcet funkˇcn´ıch hodnot.

(21)

F Pˇr´ıklad

Runge-Kuttovou metodou 4. ˇr ádu ˇreˇste poˇc áteˇcn´ı úlohu y⁰=t²− y , y (0) = 1, s krokem h = 0.1 na intervalu h0; 0.5i.

Re ˇsen´ıˇ Zn ´ame t0=0, y0=1, f (t, y ) = t²− y , budeme poˇc´ıtat y1, t.j.

pˇribliˇznou hodnotu ˇreˇsen´ı v bod ˇe t1=0.1.

k1 = f (0; 1) = 0²− 1 k2 = f (0 +1

20.1; 1 +1

20.1(−1)) = f (0.05; 0.95) = −0.9475 k3 = f (0 +1

20.1; 1 +1

20.1(−0.9475)) = f (0.05; 0.952625) = −0.950125 k4 = f (0 + 0.1; 1 + 0.1(−0.950125)) = f (0.1; 0.9049875) = −0.8949875 y1 = y0+1

60.1(k1+2k2+2k3+k4) .

=0.9051627.

Pro srovn ´an´ı: pˇresn ´e ˇreˇsen´ı je y = −e^−t +t²− 2t + 2 a jeho hodnota y (0.1) = 0.9051626 .

Vypoˇct ˇete pˇribliˇznou a pˇresnou hodnotu ˇreˇsen´ı v bodech 0.2; 0.3; 0.4 a 0.5.

(22)

Chyba v ´ypo ˇctu

Chyba v´ypoˇctu

a) y (tn) . . . pˇresn ´e ˇreˇsen´ı v ˇcase tn, yn . . . pˇribliˇzn ´e ˇreˇsen´ı v ˇcase tn

en=yn− y (tn) . . . glob ´aln´ı chyba aproximacenebo glob ´aln´ı diskretizaˇcn´ı chyba

b) Poˇc´ıt ´ame v koneˇcn ´e aritmetice:

rn=eyn− yn. . .zaokrouhlovac´ı chyba

Chceme poˇc´ıtat yn, ale ve skuteˇcnosti poˇc´ıt ´ame tzv. numerickou aproximacieyn.

Pro Eulerovu metodu:

(1) . . . glob áln´ı chyba aproximace en, enje pˇr´ımo úm ˇern á prvn´ı mocnin ˇe h (2) . . . zaokrouhlovac´ı chyba rn

rnje nepˇr´ımo úm ˇern á prvn´ı mocnin ˇe h (3) . . . celkov á chyba pro Eulerovu metodu O zaokrouhlovac´ıch chyb ách se m ˚uˇzeme pˇresv ˇedˇcit jen opakovaným výpoˇctem s r ˚uznou pˇresnost´ı (double precission, . . . )

(23)

Richardsonova extrapolace

Reˇsme poˇc áteˇcn´ı úlohuˇ numerickou metodou ˇr ádu p. Oznaˇc´ıme y (x ) pˇresn é ˇreˇsen´ı naˇs´ı úlohy. Volme dv ˇe r ˚uzn é velikosti kroku h: h = h1a h = h2a oznaˇcme y1=y (x , h1)pˇribliˇznou hodnotu ˇreˇsen´ı v bod ˇe x s krokem h1, y2=y (x , h2)pˇribliˇznou hodnotu ˇreˇsen´ı v bod ˇe x s krokem h2. Pak

y (x ) .

= y1(x ) + C · h^p₁ (4)

y (x ) .

= y2(x ) + C · h^p₂, (5)

kde C je konstanta v obou pˇr´ıpadech stejn á, nez ávisl á na h.

Od rovnice (4) odeˇcteme rovnici (5) a dostaneme

0 = y2− y1+C · h^p₂− C · h^p₁ ⇒ C = y2− y1

h^p₁− h^p₂. Konstantu C dosad´ıme do rovnice (5):

y (x ) .

=y2(x ) + y2− y1

h^p₁− h^p₂ · h^p₂ ⇒ y (x ) .

=y12= y2

h₁ h₂

p

− y1

_h

1 h₂

p

− 1 .

y12. . . Richardsonova extrapolace ˇreˇsen´ı y z´ıskan ´a z hodnot y1a y2.

(24)

Line ´arn´ı v´ıcekrokov ´e metody

Line ´arn´ı v´ıcekrokov ´e metody

Jednokrokov é metody . . . k urˇcen´ı yn+1potˇrebujeme jen informaci z pˇredchoz´ı ˇcasov é úrovn ˇe yn

V´ıcekrokov é metody . . . k vyˇc´ıslen´ı yn+1potˇrebujeme zn át v´ıce ˇcasových

úrovn´ı (napˇr. nestaˇc´ı zaˇc´ıt z poˇc áteˇcn´ı podm´ınky) Uvaˇzujme tˇri po sob ˇe jdouc´ı ˇcasov é úrovn ˇe

tn−1, tn=tn−1+h , tn+1=tn−1+2h

a integrujme diferenci ´aln´ı rovnici od tn−1do tn+1pomoc´ı Simpsonova pravidla:

Zb a

f (x )dx ≈h

3(f (x0) +4f (x1) +2f (x2) +4f (x3) +2f (x4) + · · · +4f (xn−1) +f (xn)) Uv ˇedomme si tak ´e, ˇze

Z t_n+1 t_n−1

y⁰(t)dt = y (tn+1) −y (tn−1) .

(25)

Tedy

y (tn+1) =y (tn−1) + Z t_n+1

t_n−1

f (t, y (t))dt ≈

≈ y (tn−1) +1

3h (f (tn−1,y (tn−1)) +4f (tn,y (tn)) +f (tn+1,y (tn+1))) Necht’ yn .

=y (tn). Dostaneme metodu

yn+1=yn−1+1

3h (f (tn−1,yn−1) +4f (tn,yn) +f (tn+1,yn+1)) .

(26)

Necht’ je nyn´ı d ´ano rovnom ˇern ´e d ˇelen´ı s krokem h:

tn,tn+1=tn+h, tn+2=tn+2h, . . . Obecn á line árn´ı k −krokov á metoda:

αkyn+k+ αk −1yn+k −1+ · · · + α0yn=h(βkfn+k+ βk −1fn+k −1+ · · · + β0fn) , kde konstanty αj, βj∈ R, α^k 6= 0 a α²0+ β₀²>0; fn .

=f (tn,y (tn)).

Je-liβk =0 =⇒ yn+klze vypoˇc´ıtat explicitn ˇe ze znalosti yn, . . . ,yn+k −1a z vyˇc´ıslen´ı funkce v pˇredchoz´ıch ˇcasov´ych ´urovn´ıch =⇒

explicitn´ı k −krokov ´a metoda

Je-liβk 6= 0, pak se yn+k vyskytuje na obou stran ´ach rovnice, a tedy metoda jeimplicitn´ı .

Pozn ámka line árn´ı– ve formuli se vyskytuj´ı jen line árn´ı kombinace yna f (tn,yn).

(27)

Pˇr´ıklady

Ctyˇrkrokov ´a line ´arn´ıˇ explicitn´ı Adams–Bashforthova metoda yn+4=yn+3+ 1

24h (55fn+3− 59fn+2+37fn+1− 9fn) Ctyˇrkrokov ´a line ´arn´ıˇ implicitn´ı Adams–Moultonova metoda

yn+4=yn+3+ 1

24h (9fn+4+19fn+3− 5fn+2− 9fn+1) Pozn ámka Neˇz m ˚uˇzeme aplikovat k −krokovou metodu, mus´ıme zn át k poˇc áteˇcn´ıch hodnot y0, . . . ,yk −1, kde y0je dan á poˇc áteˇcn´ı podm´ınka, y1, . . . ,yk −1mus´ı být n ˇejak vypoˇcteny, napˇr. Eulerovou metodou nebo RK metodou. V kaˇzd ém pˇr´ıpad ˇe budou data obsahovat numerick é chyby a je d ˚uleˇzit é v ˇed ˇet, jak to ovlivn´ı dalˇs´ı aproximace yn, n ≥ k , kter é poˇc´ıt áme k −krokovou metodou. Tedy zaj´ım á n ásstabilitanumerick é metody vzhledem k malým perturbac´ım poˇc áteˇcn´ıch dat.

(28)

F Stabilita k−krokov ´ych metod

Jak zjistitstabilitu ?

{tn} . . . rovnom ˇern é d ˇelen´ı s krokem h . Obecn á line árn´ı k −krokov á metoda:

α0yn+α1yn+1+· · ·+αkyn+k =h (β0f (tn,yn) + β1f (tn+1,yn+1) + · · · + βkf (tn+k,yn+k)) , kde α0, . . . , αk a β0, . . . , βk jsou re ´aln ´e konstanty, αk 6= 0, α²0+ β₀²>0 .

Oznaˇcme polynomy

ρ(z) =

k

X

j=0

αjz^j = α0+ α1z + · · · + αkz^k 1. charakteristick´y polynom

σ(z) =

k

X

j=0

βjz^j= β0+ β1z + · · · + βkz^k 2. charakteristick´y polynom

(29)

F

V ˇetaPodm´ınka stability metody

Line árn´ı v´ıcekrokov á metoda je numericky stabiln´ı pro libovolnou diferenci áln´ı rovnici y⁰=f (t, y ), kde f je Lipschitzovsk á,

pr ´av ˇe kdyˇz

koˇreny 1. charakteristick ého polynomu ρ(z) leˇz´ı uvnitˇr jednotkov ého uzavˇren ého kruhu, pˇriˇcemˇz koˇreny leˇz´ıc´ı na jednotkov é kruˇznici jsou jednon ásobn é.

Pozn ´amka f je Lipschitzovsk ´ana oblasti J × D ⇐⇒

∃ L > 0 : |f (t, y ) − f (t, z)| ≤ L|y − z| ∀(t, y ), (t, z) ∈ J × D .

Pozn ´amka Nezkoumali jsme chybu aproximace, t.j. pˇresnost k −krokov´ych metod.

(30)

F

Pˇr´ıklady

1. Adams–Bashforthova metoda yn+4=yn+3+ 1

24h (55fn+3− 59fn+2+37fn+1− 9fn) yn+4− yn+3= 1

24h (55fn+3− 59fn+2+37fn+1− 9fn) Tedy pro 1. charakteristicky polynom ρ(z) m ´ame

ρ(z) = z⁴− z³=z³(z − 1) = 0 =⇒

z = 0 je trojn ásobný koˇren uvnitˇr jednotkov ého kruhu z = 1 leˇz´ı na jednotkov é kruˇznici, je jednon ásobný

=⇒ metoda je numericky stabiln´ı . 2. Tˇr´ıkrokov ´a metodaˇr ´adu 6

11yn+3+27yn+2− 27yn+1− 11yn=3h (fn+3+9fn+2+9fn+1+fn) ρ(z) = 11z³+27z²− 27z − 11 = 0 (reciprok ´a rovnice) Koˇreny: z1=1, z2 .

= −0, 3189, z3 .

= −3, 1356 =⇒ |z3| > 1 ,

=⇒ metoda nen´ı numericky stabiln´ı.

(31)

F

3. Urˇcete vˇsechny hodnoty b ∈ R, pro kter é je line árn´ı k−krokov á metoda yn+3+ (2b − 3)(yn+2− yn+1) −yn=hb(fn+2+fn+1)

numericky stabiln´ı.

ρ(z) = z³+ (2b − 3)(z²− z) − 1 = 0 Protoˇze ρ(1) = 0 je jeden koˇren z = 1, a to jednon ´asobn´y,

z³+ (2b − 3)(z²− z) − 1 = (z − 1) · (z²+z + 1 + z(2b − 3)

| {z }

=0 := ρ1(z) = z²+z(2b − 2) + 1 Zbýv á naj´ıt koˇreny ρ1(z). Vyzkouˇs´ıme ρ1(1) = 2b =⇒ b 6= 0 , jinak by z = 1 nebyl jednon ásobný a metoda by nebyla stabiln´ı. Ze stejn ého d ˚uvodu, protoˇze ρ(−1) = −2b + 4, mus´ı být b 6= 2.

(32)

F

Jak ´e jsou tedy koˇreny ρ1? Oznaˇcme je z1,z2. Pak

(z − z1)(z − z2) =z²+z(2b − 2) + 1 =⇒ −(z1+z2)z +z1z2=z(2b − 2)+1. Tedy z1z2=1 . Ale z16= ±1 , z26= ±1, tedy z1,z2jsou komplexn´ı.

D = 4(b − 1)²− 4 < 0 ⇐⇒ b ∈ (0, 2) . Z ´av ˇer

Je-li b ∈ (0, 2), jsou koˇreny ρ(z):

z1=1, z2,3=1 − b ±ip

1 − (b − 1)², z26= z3, |z2,3| < 1 , a tedy vˇsechny koˇreny ρ(z) leˇz´ı pro b ∈ (0, 2) v uzavˇren ´em jednotkov ´em

kruhu =⇒

metoda je numericky stabiln´ı ⇐⇒ b ∈ (0, 2) .

(33)

F

V ˇeta Nutn ´a podm´ınka (nikoliv postaˇcuj´ıc´ı) pro konvergenci v´ıcekrokov ´e metody je, ˇze metoda je numericky stabiln´ı.

Line árn´ı k −krokov á metoda, jej´ıˇz charakteristický polynom je ρ(z) = z^k− z^{k −1}. . . tzv.Adamsovy metody

– explicitn´ı . . . Adamsovy–Bashforthovy metody – implicitn´ı . . . Adamsovy–Moultonovy metody

Line árn´ı k −krokov á metoda, jej´ıˇz charakteristický polynom je ρ(z) = z^k− z^{k −2}:

– explicitn´ı . . . Nystr ¨omova metoda

– implicitn´ı . . . Milneova–Simpsonova metoda.

Pozn ámka Lze rovn ˇeˇz zkoumat tzv.absolutn´ı stabilitu (A–stabilitu)line árn´ıch v´ıcekrokových metod. A–stabilitou se nebudeme zabývat.

(34)

”Stiff”syst ´emy

”Stiff”rovniceje diferenci áln´ı rovnice, pro kterounumerick á metoda je numericky nestabiln´ı, pokud nen´ı krok metody extr émn ˇe malý.

V rovnici se vyskytuj´ı ˇcleny, kter ´e zp ˚usobuj´ı rychlou zm ˇenu ˇreˇsen´ı. Jsou to napˇr´ıklad rovnice typu

y⁰=ky + f (t), kde k ∈ C, |k| velk ´a nebo soustavy

y⁰=Ky + f(t) ,

kde K m á jedno vlastn´ı ˇc´ıslo λ ∈ C takov é, ˇze |λ| je velk á v porovn án´ı s f(t) nebo <λi <0, 1 ≤ i ≤ n, ale

max

1≤i≤n|<λi| >> min

1≤i≤n|<λi| . M´ıra stiff

R = max |<λi|

min |<λi| , λi je vlastn´ı ˇc´ıslo Jacobivy matice dan ´eho syst ´emu.

Zat´ım ˇz ádn á, obecn ˇe akceptovan á definice ”stiff”neexistuje.

(35)

Metody typu prediktor–korektor

Oznaˇcme

AB – Adamsova–Bashforthova metoda (explicitn´ı k −krokov á 2. ˇr ádu) AM – Adamsova–Moultonova metoda (implicitn´ı k −krokov á 2. ˇr ádu) Idea:

Prediktor– v naˇsem pˇr´ıpad ˇe explicitn´ı metoda (AB) . Jej´ı v´ystup povaˇzujeme zameziv´ysledek

˜

yn+2=yn+1+h

2(3f (tn+1,yn+1) −f (tn,yn)) .

Nyn´ı”oprav´ıme”(correct)tuto hodnotu pomoc´ı implicitn´ı metody (AM), pˇriˇcemˇz meziv´ysledek ˜yn+2pouˇzijeme na prav ´e stran ˇe. Dostaneme

yn+2=yn+1+h

2(f (tn+1,yn+1) +f (tn+2, ˜yn+2)) .

(36)

Literatura k dal ˇs´ımu studiu

Fajmon B., R ˚uˇziˇckov ´a I.: Matematika 3. Fakulta elektrotechniky a komunikaˇcn´ıch technologi´ı, VUT Brno, 2011.

Kub´ıˇcek M., Dubcov á M., Janovsk á D.: Numerick é metody a algoritmy, V ˇSCHT Praha, 2005 (second edition).

Rasmuson A., Andersson B., Olsson L., Andersson R.: Mathematical Modeling in Chemical Engineering. Cambridge University Press, 2014.

Recktenwald G.: Numerical Integration of Ordinary Differential Equations for Initial Value Problems. Portland State University 2007,

gerry@me.pdx.edu.

Turz´ık D. et all: Matematika II, V ˇSCHT Praha, 2002.