Numerická matematika. Úvodní informace. Viz Kontakt: Petr Sváček, KN:D 201

Numerick´a matematika ´ Uvodn´ ı informace

Viz http://mat.fs.cvut.cz I

Kontakt: Petr Sv´aˇcek, KN:D 201

Soustavy rovnic, souvislost s prax´ı

I

Tˇeleso nahrad´ıme diskr´etn´ımi body, hled´ame nezn´am´e fyzik´aln´ı veliˇciny v tˇechto bodech.

I

Uˇzijeme fyzik´aln´ı vztahy, z´ısk´ame soustavu neline´arn´ıch/line´arn´ıch rovnic.

I

Hled´ame ˇreˇsen´ı t´eto soustavy: n - poˇcet nezn´am´ych.

Soustavy line´arn´ıch rovnic     

a11 a12 . . . a1n a21 a22 . . . a2n .. .. . . an1 an2 . . . ann

      ·  

x1 x2 .. .





    =  

xn

b1 b2 .. .

    

bn

I

Zn´ame: transponovan´a/symetrick´a, hodnost, regul´arn´ı/singul´arn´ı, Frobeniova vˇeta.

I

Zn´ame: Gaussova eliminace, LU rozklad, determinant, inverzn´ı matice.

I

Zn´ame: Vlastn´ı ˇc´ıslo a vektor. Au = λu,

I

u 6= 0,

Nov´ e: Velk´ e ˇr´ıdk´ e matice, efektivita ˇreˇsen´ı.

Pˇr´ım´e metody ˇreˇsen´ı

Ax = b Pˇr´ım´e metody jsou takov´e metody, kter´e vedou k nalezen´ı pˇresn´eho ˇreˇsen´ı v pˇredem zn´am´em poˇctu krok˚ u. V´yhody: Vˇzdy vedou k v´ysledku. Nev´yhody: Rychle rostouc´ı v´ypoˇcetn´ı a pamˇet’ov´a n´aroˇcnost operac´ı ≈ n3 , pamˇet’ ≈ n2 . I

Gaussova eliminace.

I

Inverzn´ı matice.

I

Cramerovo pravidlo.

I

LU rozklad.

LU rozklad Postup ˇreˇsen´ı

I

Uvaˇzujme soustavu s matic´ı A = LU typu n × n LUx = b

I

ˇ s´ıme soustavu ≈ n2 . Reˇ Ly = b

I

ˇ s´ıme soustavu ≈ n2 . Reˇ Ux = y

I

Jak z´ıskat LU rozklad? (2 moˇznosti, viz cviˇcen´ı)

Shrnut´ı: Pˇr´ım´e metody ˇreˇsen´ı

Pˇr´ım´e metody ˇreˇsen´ı Ax = b I

Zaruˇceno, ˇze dos´ahnou v´ysledku.

I

Rostouc´ı v´ypoˇcetn´ı resp. pamˇet’ov´a n´aroˇcnost s n3 resp. n2 .

Efektivnˇejˇs´ı ˇreˇsen´ı? I

Staˇc´ı pˇribliˇzn´e ˇreˇsen´ı, hled´ame vhodn´y pˇredpis X (k+1) = F (X (k) , . . . )

I

viz posloupnost xn+1 =

1 2 ( + xn ) 2 xn

Pˇr´ıklad Soustavu Ax = b je moˇzn´e pˇrepsat jako x = Ux + V , Zvol´ıme X (0) , spoˇc´ıt´ame X (0) , X (1) , X (2) , ... Pokud konverguje nez´avisle na X (0) , ˇr´ık´ame, ˇze odpov´ıdaj´ıc´ı metoda je konvergentn´ı. Jak mˇeˇrit vzd´alenost X (k) od X ∗ ? Norma. Je-li V = 0, pak X (n) z´avis´ı na velikosti U (normˇe) X (n) = Un X (0)

Prost´a iteraˇcn´ı metoda

Je-li soustava tvaru x = Ux + V , Iteraˇcn´ı metoda X (k+1) = UX (k) + V , Pokud konverguje, konverguje k ˇreˇsen´ı. Konvergence z´avis´ı na vlastnostech matice U, pˇr´ıpadnˇe na vlastnostech matice A, z kter´e vznikla.

Velikost vektoru. Norma  

 3 x =  −4  , 0

  e=  

0.03 −0.04 0.03 −0.01 −0.01

     

Norma kxk =?, kek =? Norma vektoru k · k : Rn → R, α ∈ R kαuk ≤ |α|kuk ku + v k ≤ kuk + kv k kuk = 0 ⇔ u = 0 Norma matice lze definovat pomoc´ı normy vektoru kAxk x6=0 kxk

kAk = sup kAxk = sup x,kxk=1

Norma matice a vektoru I

ˇ adkov´a norma R´ kAkm = max i

I

|ai,j |

k~x km = max |xi |

|ai,j |

k~x k` =

j

m X i=1

I

|xi |

k~x kE =

X n

|xi |2

i=1

Spektrum a spektr´aln´ı polomˇer. Pozn. spektr´aln´ı norma matice kAk2 . Vztah normy matice a vektoru kAxk ≤ kAkkxk

I

n X i=1

Frobeniova/Euklidovsk´a norma X 1 m X n 2 kAkF = |ai,j |2 i=1 j=1

I

i

j=1

Sloupcov´a norma kAk` = max

I

n X

Vztah normy a spektr´aln´ıho polomˇeru.

1

2

Vlastnosti matice Matice A je I diagon´ alnˇ e dominantn´ı (DD), pokud ∀i X |aii | ≥ |aij | j6=i I

ostˇre diagon´ alnˇ e dominantn´ı (ODD), pokud ∀i X |aii | > |aij | j6=i

I

pozitivnˇ e semidefinitn´ı pokud pro libovoln´y vektor plat´ı x T Ax =

n X

xi aij xj ≥ 0

i,j=1 I I I

pozitivnˇ e definitn´ı - ostr´a nerovnost pro nenulov´y vektor symetrick´ a pozitivnˇ e definitn´ı (SPD) Symetrick´e pozitivnˇe definitn´ı matice ⇔ kladn´a vlastn´ı ˇc´ısla nebo kladn´e hlavn´ı minory.

Vlastnosti matice Matice A je I

reducibiln´ı/rozloˇ ziteln´ a, pokud ∃ P permutaˇcn´ı matice tak, ˇze   B C T P AP = 0 D

I

ireducibiln´ı/nerozloˇ ziteln´ a, pokud nen´ı rozloˇziteln´a/reducibiln´ı.

I

Graf matice: uzly 1, . . . , n, orientovan´e hrany od i k j pro aij 6= 0

I

Graf - silnˇ e souvisl´ y, od libovoln´eho i se dostanu k libovoln´emu j. Nerozloˇ ziteln´ a matice m´ a silnˇ e souvisl´ y graf.

I

je ireducibilnˇ e/nerozloˇ zitelnˇ e diagon´ alnˇ e dominantn´ı (IDD), pokud (i) je diagon´alnˇe dominantn´ı, (ii) ∃i X |aii | > |aij | j6=i

(iii) je ireducibiln´ı, I

Souvislost s regularitou matice.

Prost´a iteraˇcn´ı metoda Je-li soustava tvaru x = Ux + V , Iteraˇcn´ı metoda X (k+1) = UX (k) + V , Plat´ı I

Nutn´ a a postaˇ cuj´ıc´ı podm´ınka: Prost´a iteraˇcn´ı metoda je konvergentn´ı pr´avˇe tehdy, kdyˇz ρ(U) < 1.

I

Postaˇ cuj´ıc´ı podm´ınka: kk - maticova norma. Je-li kUk < 1, pak prost´a iteraˇcn´ı metoda je konvergentn´ı.

I

Odhad chyby: kx (k) − x ∗ k ≤

kUk kx (k) − x (k−1) k, 1 − kUk

kx (k) − x ∗ k ≤

kUkk kx (1) − x (0) k, 1 − kUk

Prost´a iteraˇcn´ı metoda: Pˇr´ıklad. Pˇr. 1. Urˇcete vˇsechna p ∈ R, pro kter´a konverguje prost´a iteraˇcn´ı metoda pro soustavu tvaru X = UX + V , kde V = (−1, 2, 0.5)T ,   −0.8 p 2 0 1 0.8 0  , U= 1 −1 p Pro p = 0.2, X (0) = V urˇcete X (1) touto metodou. Pˇr. 2. Je d´ana soustava rovnic ~x = R~x + ~s , kde ~s = (1.7, −1.8, 0.5)T   0.1 t + 1 −0.3 0.1  R =  0.4 −0.3 0.6 0.1 t + 0.5 Urˇcete vˇsechna t ∈ R, pro nˇeˇz je splnˇena postaˇcuj´ıc´ı podm´ınka konvergence prost´e iteraˇcn´ı metody pro danou soustavu. Pro t = −0.5 urˇcete ~x (2) touto metodou pˇri volbˇe ~x (0) = ~0. Vypoˇctˇete k~x (2) − ~x (1) km . Odhadnˇete chybu.

Pˇredn´aˇska ˇc. 2

I

Jacobiho iteraˇcn´ı metoda, postup v´ypoˇctu, krit´eria konvergence.

I

Gauss-Seidelova iteraˇcn´ı metoda, postup v´ypoˇctu, krit´eria konvergence.

I

Super-relaxaˇcn´ı metoda (SOR).

Opakov´an´ı

I

Norma matice - ˇr´adkov´a, sloupcov´a, Euklidovsk´a.

I

PD, SPD, DD, nerozloˇzitelnˇe/ireducibilnˇe diagon´alnˇe dominantn´ı (IDD), ODD (ˇr´adky nebo sloupce), regularita matice

I

Prost´a iteraˇcn´ı metoda.

I

Souvislost pojm˚ u

I

Dopln´ıme ˇc´ıslo podm´ınˇennosti.

ˇ ıslo podm´ınˇenosti a ˇspatnˇe podm´ınˇenn´a matice C´ Doplnˇen´ı pojm˚ u I I I

Je-li det A → 0, nemus´ı b´yt soustava jeˇstˇe ˇspatnˇe ˇreˇsiteln´a. Vad´ı tzv. ˇspatnˇe podm´ınˇen´a matice (tj. velk´e κ(A)) ˇ ıslo podm´ınˇ C´ enosti κ(A) Vezmeme ˇreˇsen´ı x ∗ a x˜ pro bl´ızkou ˜ pravou stranu b, Ax ∗ = b

A˜ x = b˜

Uk´aˇzeme, ˇze plat´ı ˜ k(x ∗ − x˜)k kb − bk −1 ≤ kA kkAk | {z } kbk kx ∗ k κ(A)

Pro SPD matice plat´ı κ2 (A) = λMAX /λMIN Pˇr.   1 1 A= , b = (2, 1.99999)T , b˜ = (2, 1.999999)T . 1 0.99999

Prost´a iteraˇcn´ı metoda Soustava tvaru x = Ux + V, Iteraˇcn´ı metoda X(k+1) = UX(k) + V, Plat´ı I

Nutn´ a a postaˇ cuj´ıc´ı podm´ınka: Prost´a iteraˇcn´ı metoda je konvergentn´ı pr´avˇe tehdy, kdyˇz ρ(U) < 1.

I

Postaˇ cuj´ıc´ı podm´ınka: Je-li kUk < 1, pak prost´a iteraˇcn´ı metoda je konvergentn´ı. Nav´ıc plat´ı kX(k) − x∗ k ≤

kUk kX(k) − X(k−1) k, 1 − kUk

kX(k) − x∗ k ≤

kUkk kX(1) − X(0) k, 1 − kUk

Prost´a iteraˇcn´ı metoda: Pˇr´ıklad. Pˇr. 1. Urˇcete vˇsechna p ∈ R, pro kter´a konverguje prost´a iteraˇcn´ı metoda pro soustavu tvaru X = UX + V , kde V = (−1, 2, 0.5)T ,   −0.8 p 2 0 1 0.8 0  , U= 1 −1 p Pro p = 0.2, X (0) = V urˇcete X (1) touto metodou. Pˇr. 2. Je d´ana soustava rovnic ~x = R~x + ~s , kde ~s = (1.7, −1.8, 0.5)T   0.1 t + 1 −0.3 0.1  R =  0.4 −0.3 0.6 0.1 t + 0.5 Urˇcete vˇsechna t ∈ R, pro nˇeˇz je splnˇena postaˇcuj´ıc´ı podm´ınka konvergence prost´e iteraˇcn´ı metody pro danou soustavu. Pro t = −0.5 urˇcete ~x (2) touto metodou pˇri volbˇe ~x (0) = ~0. Vypoˇctˇete k~x (2) − ~x (1) km . Odhadnˇete chybu.

Jak z matice A z´ısk´am matici U? Odvozen´ı klasick´ych iteraˇcn´ıch metod

Vezmeme postupnˇe rovnici i = 1, 2, 3, . . . , n a z t´e nalezneme xi : X X aij xj + aii xi + aij xj = bi j
j>i

Vyj´adˇr´ıme xi  xi =

1  bi − aii

 X j
aij xj −

X

aij xj 

j>i

Pozn. lze t´eˇz v maticov´em z´apise A = L + D + P, dopln´ıme iterace (k+1) (k/k+1) xi , xj (Jacobi, Gauss-Seidel).

Jacobiho iteraˇcn´ı metoda Pro A = L + D + P. Dost´av´ame Dx(k+1) = −(L + P)x(k) + b, a Jacobiho iteraˇcn´ı metoda −1 x(k+1) = −D−1 (L + P) x(k) + D | {z b} . {z } | UJ

VJ

Plat´ı I

Je-li matice A ODD(IDD), pak Jacobiho iteraˇcn´ı metoda je konvergentn´ı.

I

Jacobiho iteraˇcn´ı metoda je konvergentn´ı pr´avˇe tehdy, kdyˇz ρ(UJ ) < 1.

I

Vlastn´ı ˇc´ısla matice UJ jsou koˇreny rovnice det(L + λD + P) = 0

Jacobiho iteraˇcn´ı metoda: Pˇr´ıklad. Pˇr. 3. Je d´ana soustava  1  −1 A= 0

Ax = b, kde  p 0 2 −1  , p 4



 1 b= 0  4

a) Pro jak´e hodnoty parametru p ∈ R je matice A ODD (IDD, SPD)? Zd˚ uvodnˇete! b) Urˇcete vˇsechny hodnoty parametru p ∈ R, pro kter´e je Jacobiho iteraˇcn´ı metoda pro danou soustavu konvergentn´ı. c) Volte p = −1 a X (0) = (1, −2, 1)T . Spoˇctˇete Jacobiho iteraˇcn´ı metodou X (1) .

Kde se pouˇz´ıvali/j´ı iteraˇcn´ı metody?

V´ystavba pˇrehrady Orl´ık, 1954-1961. Betonov´a hr´az 450 m dlouh´a, v´yˇska 91 m.

Sloˇzkov´y z´apis klasick´ych iteraˇcn´ıch metod

I

Jacobiho iteraˇ cn´ı metoda   X X 1  (k) (k) bi − aij xj − aij xj  xi (k+1) = aii j
I

j>i

Gauss-Seidelova iteraˇ cn´ı metoda   X X 1  (k+1) (k) xi (k+1) = bi − aij xj − aij xj  aii j
j>i

Gauss-Seidelova iteraˇcn´ı metoda Pro Ax = b zap´ıˇseme A = L + D + P. Dost´av´ame (D + L)x(k+1) = −Px(k) + b, a tedy iteraˇcn´ı metodu x(k+1) = −(D + L)−1 P x(k) + (D + L)−1 b . | {z } | {z } UGS

VGS

Plat´ı I Je-li matice A ODD (IDD), pak je GS iteraˇ cn´ı metoda konvergentn´ı. I Je-li matice A SPD, pak je GS iteraˇ cn´ı metoda konvergentn´ı. I GS iteraˇ cn´ı metoda je konvergentn´ı pr´avˇe tehdy, kdyˇz ρ(UGS ) < 1. I Vlastn´ ı ˇc´ısla matice UGS jsou koˇreny rovnice det(λL + λD + P) = 0

Iteraˇcn´ı metody: Pˇr´ıklad. Pˇr. 4. D´ana soustava line´arn´ıch rovnic AX = B, kde     1 2 0 52 A =  0 4 1 , B =  p , −2 p 1 5 a) Urˇcete vˇsechny hodnoty parametru p ∈ R, pro kter´e je splnˇena nˇekter´a z postaˇcuj´ıc´ıch podm´ınek konvergence Gauss-Seidelovy iteraˇcn´ı metody. b) Urˇcete vˇsechny hodnoty parametru p ∈ R, pro nˇeˇz je splnˇena nutn´a a postaˇcuj´ıc´ı podm´ınka Gauss-Seidelovy iteraˇcn´ı metody. c) Volte p = 0 a spoˇctˇete X (1) pokud X (0) = B. Pˇr. 5. D´ana soustava 

5  1 2

1 2 0

   p 2 0  · X =  p , 1 −2

a) Urˇcete vˇsechna p ∈ R, pro nˇeˇz je splnˇena nˇekter´a z postaˇcuj´ıc´ıch podm´ınek (pro matici A, ne UG ) GS iteraˇcn´ı metody. b) Urˇcete vˇsechna p ∈ R, pro nˇeˇz je splnˇena nutn´a a postaˇcuj´ıc´ı podm´ınka GS iteraˇcn´ı metody. c) Volte p = 0 a spoˇctˇete X (1) pokud X (0) = B.

Superrelaxaˇcn´ı metoda (SOR) I

SOR (ω - relaxaˇcn´ı parameter) xi (k+1) =

(k) xi

  X X ω  (k+1) (k)  + aij xj − aij xj bi − aii j
j≥i

nebo xi (k+1)

  X X ω  (k+1) (k) (k) = bi − aij xj − aij xj  + (1 − ω)xi aii j
I I

Dxk+1 = Dxk + ω(b − Lxk+1 − Dxk − Uxk ) Plat´ı I I I

I

j>i

Je-li matice A ODD, ω ∈ (0, 1i, pak je SOR konvergentn´ı. Je-li matice A SPD, ω ∈ (0, 2), pak je SOR konvergentn´ı. Je-li SOR konvergentn´ı, pak ω ∈ (0, 2).

ρ(USOR ) je spojit´a funkce v ω.

Superrelaxaˇcn´ı metoda (SOR)

2 -1 0 0 0 0 0 0

-1 2 -1 0 0 0 0 0

0 -1 2 -1 0 0 0 0

0 0 -1 2 -1 0 0 0

0 0 0 -1 2 -1 0 0

0 0 0 0 -1 2 -1 0

0 0 0 0 0 -1 2 -1

0 0 0 0 0 0 -1 2

Optim´aln´ı parameter (za urˇcit´ych pˇredpoklad˚ u, ρ = ρ(UGS )) ωopt =

2 p 1 + 1 − ρ2

Pro dan´y obr´azek: 20 iterac´ı GS (ω = 1) 0.9220 ≈ 0.18 20 iterac´ı SOR (ω = 1.6) 0.620 ≈ 3.6 × 10−5

Pˇr´ıklad. SOR Pˇr. 6. Je d´ana soustava  1  −1 A= 0

Ax = b, kde  p 0 2 −1  , p 4



 1 b= 0  4

a) Pro jak´e hodnoty parametru p ∈ R je matice A SPD? Zd˚ uvodnˇete! b) Volte p = −1 a X (0) = (1, −2, 1)T . Volte ω = 1.5 a spoˇctˇete SOR iteraˇcn´ı metodou X (1) . c) Je pro tento pˇr´ıpad SOR konvergentn´ı? Zd˚ uvodnˇete.

Pˇredn´aˇska ˇc. 3

I

Gradientn´ı metody, metoda nejvˇetˇs´ıho sp´adu, pouˇzit´ı na soustavy s SPD matic´ı.

I

Metoda sdruˇzen´ych gradient˚ u.

I

Uvaˇzujeme kvadratickou funkci g (x1 , x2 ) = 2x12 − x1 x2 + 4x22

Gradientn´ı metoda I

Pˇredpoklady: Matice A symetrick´a a pozitivnˇe definitn´ı. Hled´ame minimum funkce: 1 F (x) = xT Ax − xT b 2

Gradientn´ı metoda I

I

Pˇredpoklady: Matice A symetrick´a a pozitivnˇe definitn´ı. Hled´ame minimum funkce: 1 F (x) = xT Ax − xT b 2 Ekvivalence: ( Ax ∗ = b

I

)

⇔

(F (x∗ ) ≤ F (x)

pro libovoln´e x) .

Iteraˇcn´ı proces: x(k+1) = x(k) + αp

I

Volba smˇeru (nejvˇ etˇs´ı sp´ ad): p = − grad F (x(k) ) = b − Ax(k) .

I

Optim´ aln´ı krok ve smˇeru p
p · (b − Ax) α = argmin F x(k) + αp = p · (Ap) α

Pˇr´ıklad. Nejvˇetˇs´ı sp´ad Pˇr. 6. Je d´ana soustava  1  −1 A= 0

Ax = b, kde  −1 0 2 −1  , −1 4



 1 b= 0  4

a) Ukaˇzte, ˇze pro danou matici je ˇreˇsen´ı soustavy rovnic ekvivalentn´ı minimalizaci funkce F (x) = 21 xT Ax − xT b b) Odvod’te podm´ınku pro volbu smˇeru p tak aby pro iteraˇcn´ı proces x(k+1) = x(k) + αp doˇslo k poklesu funkce F (x) pro nˇejak´a α > 0. Zvolte smˇer, v kter´em je pokles funkce nejvˇetˇs´ı. c) Pro dan´y smˇer odvod’te volbu optim´aln´ıho kroku α, tj. tak aby pokles v dan´em smˇeru byl nejvˇetˇs´ı moˇzn´y. d) Volte X (0) = (1, −2, 1)T a spoˇctˇete metodou nejvˇetˇs´ıho sp´adu X (1) .

Metoda sdruˇzen´ych (konjugovan´ych) gradient˚ u Jak volit l´epe smˇer?

I

ortogon´ aln´ı smˇery - nevede ke zlepˇsen´ı,

I

A-ortogon´ aln´ı: pi · (Apj ) = pj · (Api ) = δij .

I

Gram-Schmidt˚ uv OG proces(?),

I

Zkus´ım: OG jen 2 posledn´ı smˇery Vol´ıme x0 , spoˇcteme p0 = b − bAx0 D´ale poˇc´ıt´ame pro k = 0, 1, ...

I

I I I I

Krok αk optim´aln´ı ve smˇeru pk , Novou iteraci xk+1 = xk + αk pk , Reziduum rk+1 = b − Axk+1 , Nov´ y smˇ er: kolmou ˇc´ast rk+1 k pk .

Jak volit l´epe smˇer? Metoda sdruˇzen´ych (konjugovan´ych) gradient˚ u I

Pˇredpoklady: Matice A symetrick´a a pozitivnˇe definitn´ı. x

(k+1)

=x

(0)

+

k X

αj pj

j=0 I

Oznaˇcme chybu ej = xj − x∗ .

I

Optim´ aln´ı krok lze zapsat αj = −

I

Tedy: Jsou-li smˇery A-ortogon´aln´ı, pak plat´ı: pi · (Aek+1 ) = 0

I

pj · (Ae0 ) pj · (Apj )

D˚ usledek?

pro i ≤ k

Shrnut´ı - soustavy line´arn´ıch rovnic Pˇr´ıklady k zamyˇslen´ı

I

Jednoduch´ e uk´ azat: I. II. III. IV. V. VI.

ρ(A) ≤ kAk, kAk2 ≤ kAkF kUk < 1 pak PIM konvergentn´ı. Uk´azat odhad chyby. V´ypoˇcet vlastn´ıch ˇc´ısel UJ a UGS, ODD v ˇr´adc´ıch, pak Jacobi konvergentn´ı (jednoduch´e), SPD pak aii > 0 a vlastn´ı ˇc´ısla jsou re´aln´a kladn´a ODD ve sloupc´ıch, pak Jacobi konvergentn´ı,

Shrnut´ı - soustavy line´arn´ıch rovnic Pˇr´ıklady k zamyˇslen´ı

I

Jednoduch´ e uk´ azat: I. II. III. IV. V. VI.

I

ρ(A) ≤ kAk, kAk2 ≤ kAkF kUk < 1 pak PIM konvergentn´ı. Uk´azat odhad chyby. V´ypoˇcet vlastn´ıch ˇc´ısel UJ a UGS, ODD v ˇr´adc´ıch, pak Jacobi konvergentn´ı (jednoduch´e), SPD pak aii > 0 a vlastn´ı ˇc´ısla jsou re´aln´a kladn´a ODD ve sloupc´ıch, pak Jacobi konvergentn´ı,

Obt´ıˇ znˇ ejˇs´ı: I I I I I I

kρ(U)k < 1 pak PIM konvergentn´ı (Jordan˚ uv tvar), ODD ve sloupc´ıch/ˇr´adc´ıch, pak GS je konvergentn´ı, IDD pak Jacobi/GS je konvergentn´ı, (viz Gerschgorin) ODD/IDD, pak SOR konvergentn´ı pro ω ∈ (0, 1i, SPD pak GS je konvergentn´ı, SPD, pak SOR je konvergentn´ı pro ω ∈ (0, 2),

Pˇredn´aˇska ˇc. 4

I

Soustavy neline´ arn´ıch rovnic

I

Probl´emy existence a jednoznaˇcnosti ˇreˇsen´ı.

I

Iteraˇcn´ı metody: metoda prost´ e iterace, Newtonova metoda.

Soustavy neline´arn´ıch algebraick´ych rovnic

Syst´ em n-rovnic o n-nezn´ am´ ych f1 (x1 , . . . , xn ) = 0 .. . fn (x1 , . . . , xn ) = 0

I

speci´aln´ı volba → soustava line´arn´ıch rovnic

I

nen´ı zaruˇcena existence ani jednoznaˇcnost ˇreˇsen´ı

Prost´a iteraˇcn´ı metoda Soustava rovnic ve tvaru F1 (x1 , . . . , xn ) = x1 .. nebo-li x = F(x) . Fn (x1 , . . . , xn ) = xn Prost´ a iteraˇ cn´ı metoda X(k+1) = F(X(k) ), Kontraktivn´ı zobrazen´ı: F : M → M, c ∈ [0, 1) tak, ˇze kF (x) − F (y )k ≤ ckx − y k Plat´ı I Je-li F kontraktivn´ ı, pak existuje x ∗ = F (x ∗ ) (pevn´y bod, Banachova vˇeta). I 1D: |F 0 (x)| < 1, I Rn : kF 0 (x)k < 1.

Pˇr´ıklad. Prost´a iterace pro neline´arn´ı rovnice

Pˇr. Je d´ana rovnice x = F(x), kde vektorov´a funkce F m´a sloˇzky 1 2 F1 (x1 , x2 ) = ( + x1 ) 2 x1 + x2 1 3 F2 (x1 , x2 ) = ( + x2 ) 2 x1 + x2 a) Urˇcete Jacobiho matici funkce F a ovˇeˇrte, zda v bodˇe [1, 1] je jej´ı norma menˇs´ı neˇz jedna. b) Volte X (0) = (1; 1)T a stanovte aproximaci X (1) = (x (1) , y (1) ) jednoho z ˇreˇsen´ı rovnice pomoc´ı prost´e iteraˇcn´ı metody.

Soustavy neline´arn´ıch algebraick´ych rovnic Newtonova metoda pro pˇr´ıpad 1d

V bodˇe x k uˇzijeme Taylor˚ uv polynom 0 = f (x) = f(xk ) + f 0 (xk )(x − xk ) + R2 (x) Zanedb´an´ım dostaneme vzorec pro x ≈ x k+1  −1 f (x k ). x k+1 = x k − f 0 (x k ) Obecn´y vzorec pro soustavu F (x) = 0. Co je F 0 (x) ?  −1 x k+1 = x k − F 0 (x k ) F (x k ).

Newtonova metoda - odvozen´ı Odvozen´ı pro 2 neline´ arn´ı rovnice o 2 nezn´ am´ych

Oznaˇcme k−t´e pˇribl´ıˇzen´ı jako A = [x (k) , y (k) ] 0 = f (x ∗ , y ∗ ) = f (A) + fx (A) (x ∗ − x (k) ) + fy (A) (y ∗ − y (k) ) + . . . 0 = g (x ∗ , y ∗ ) = g (A) + gx (A) (x ∗ − x (k) ) + gy (A) (y ∗ − y (k) ) + . . . Zanedb´ame-li dalˇs´ı ˇcleny Taylorova rozvoje fx (A)∆x +

fy (A)∆y + f (A)

= 0,

gx (A)∆x + gy (A)∆y + g (A) = 0,

dost´av´ame soustavu line´arn´ıch rovnic pro ∆x = x (k+1) − x (k) ,

∆y = y (k+1 ) − y (k) .

Newtonova iteraˇcn´ı metoda 1. Zvol´ıme poˇc´ateˇcn´ı pˇribl´ıˇzen´ı [x 0 , y 0 ]. 2. Postupnˇe pro k = 0, 1, . . . a) Sestav´ıme soustavu rovnic v bodˇe A = [x (k) , y (k) ], fx (A)∆x + fy (A)∆y + f (A) = 0 gx (A)∆x + gy (A)∆y + g (A) = 0

b) Najdeme ˇreˇsen´ı t´eto soustavy ∆x , ∆y c) Spoˇcteme nov´e pˇribl´ıˇzen´ı  (k+1)   (k)    ∆x x x = + ∆y y (k+1) y (k) Nen´ı zaruˇ cena konvergence k ˇreˇsen´ı. Metoda m˚ uˇ ze selhat (viz b)) Konvergence z´ avis´ı na poˇ c´ ateˇ cn´ım pˇribl´ıˇ zen´ı.

Soustavy neline´arn´ıch algebraick´ych rovnic Probl´emy existence a jednoznaˇcnosti ˇreˇsen´ı

Pˇr. 9. Je d´ana soustava rovnic x2 + y2 = 1 4

y = 2 cos(πx)

a) Graficky zn´azornˇete pˇribliˇzn´y poˇcet a polohu koˇren˚ u t´eto soustavy neline´arn´ıch rovnic. b) Volte ~x (0) = (−0.5; −1)T a vypoˇctˇete ~x (1) Newtonovou metodou c) Urˇcete k~x (1) − ~x (0) k` Pˇr. 10. Je d´ana soustava rovnic 1 −y =0 2x

x 2 + 4y 2 = 4

a) Graficky zn´azornˇete pˇribliˇzn´y poˇcet a polohu koˇren˚ u t´eto soustavy neline´arn´ıch rovnic. b) Volte ~x (0) = (1; 0)T a vypoˇctˇete ~x (1) Newtonovou metodou c) Urˇcete k~x (1) − ~x (0) km

Taylor˚ uv polynom, n´ahrada derivace I

Taylor˚ uv polynom f (x + h) = Tn (h) + Rn+1 (h) =

n X f (k) (x) k=0

k!

hk + Rn+1 (h)

I

Je-li f dostateˇcnˇe hladk´a, pak |Rn+1 (h)| ≤ Chn+1 = O(hn+1 )

I

Vyuˇzijeme k n´ahradˇe derivac´ı. f (x + h) − f (x) = . . . f (x − h) − f (x) = . . . f (x + h) − f (x − h) = . . .

I

Druh´a derivace f 00 (x), pˇresnˇejˇs´ı n´ahrada f 0 (x) z hodnot v x, x − h, x − 2h.

Pˇredn´aˇska ˇc. 5 Numerick´e ˇreˇsen´ı ODR.

I

Taylor˚ uv polynom a odvozen´ı diferenˇcn´ıch n´ahrad.

I

ODR, Cauchyova u ´loha: Matematick´y popis fyzik´aln´ıch jev˚ u.

I

Obecn´a jednokrokov´a metoda. Explicitn´ı a implicitn´ı Eulerova metoda.

Obyˇcejn´a diferenci´aln´ı rovnice, Cauchyova u´loha I

Cauchyova u ´loha pro ODR 1. ˇr´adu y 0 = f (x, y ),

I

Cauchyova u ´loha pro soustavu rovnic ~y 0 = ~f (x, ~y ),

I

y (x0 ) = y0

~y (x0 ) = ~y 0 .

Cauchyova u ´loha pro rovnici n-t´eho ˇr´adu y (n) = G (x, y , y 0 , . . . , y (n−1) ), y (x0 ) = y0 , . . . y (n−1) (y0 ) = y0 ,

I

Post. podm´ınky pro existenci a jednoznaˇcnost ˇreˇsen´ı: MA3.

Taylor˚ uv polynom, n´ahrada derivace I

Taylor˚ uv polynom f (x + h) = Tn (h) + Rn+1 (h) =

n X f (k) (x) k=0

k!

hk + Rn+1 (h)

I

Je-li f dostateˇcnˇe hladk´a, pak |Rn+1 (h)| ≤ Chn+1 = O(hn+1 )

I

Vyuˇzijeme k n´ahradˇe derivac´ı. f (x + h) − f (x) = . . . f (x − h) − f (x) = . . . f (x + h) − f (x − h) = . . .

I

Druh´a derivace f 00 (x), pˇresnˇejˇs´ı n´ahrada f 0 (x) z hodnot v x, x − h, x − 2h.

Diferenˇcn´ı n´ahrady I

dopˇredn´a diference

I

f (x0 + h) − f (x0 ) = f 0 (x0 ) + O(h), h zpˇetn´a diference

I

I

I

f (x0 ) − f (x0 − h) = f 0 (x0 ) + O(h), h 1. centr´aln´ı diference f (x0 + h) − f (x0 − h) = f 0 (x0 ) + O(h2 ), 2h 2. centr´aln´ı diference f (x0 + h) − 2f (x) + f (x0 − h) = f 00 (x0 ) + O(h2 ). h2 2. zpˇetn´a diference (2. ˇr´adu), A, B, C dopoˇc´ıt´ame Af (x0 ) + Bf (x0 − h) + Cf (x0 − 2h) = f 0 (x0 ) + O(h2 ), 2h

Princip numerick´eho ˇreˇsen´ı y

3

[x3 , Y ]

[x3 , y(x3)]

y=y(x)

x x

0

x

1

x

2

x

3

x

4

I

Pˇresn´ e ˇreˇsen´ı: y = y (x) pro x ∈ ha, bi.

I

Pˇribliˇ zn´ e ˇreˇsen´ı Y i : krok h = (b − a)/n, xi = a + ih,

I

Y i ≈ y (xi )

Uˇzijeme diferenˇcn´ı n´ahrady: Eulerova metoda (explicitn´ı resp. implicitn´ı) y 0 (xn ) ≈

Y n+1 − Y n , h

y 0 (xn+1 ) ≈

Y n+1 − Y n h

Explicitn´ı a implicitn´ı Eulerova metoda. I

Explicitn´ı Eulerova metoda Y n+1 − Y n = f (xn , Y n ), h

I

Implicitn´ı Eulerova metoda Y n+1 − Y n = f (xn+1 , Y n+1 ), h

I

Y n+1 = Y n + hf (xn , Y n ),

Y n+1 −hf (xn+1 , Y n+1 ) = Y n

ˇ s´ıme u Stabilita: Reˇ ´lohu (λ ∈ C) y 0 = λy ,

y (0) = 1,

kde zn´ame analytick´e ˇreˇsen´ı. Jak se chov´a numerick´e ˇreˇsen´ı? I

Glob´ aln´ı chyba metody: O(h).

Eulerova explicitn´ı metoda - srovn´an´ı s pˇresn´ym ˇreˇsen´ım Rovnice x¨ + x = 0

I

Konvergence k ˇreˇsen´ı pro h = 2π/20, 2π/40, , 2π/80, 2π/160.

Eulerova explicitn´ı metoda - chyba v z´avislosti na h Rovnice x¨ + x = 0

I

Konvergence k ˇreˇsen´ı pro h = 2π/20, 2π/40, , 2π/80, 2π/160.

Zobecnˇen´ı Eulerovy explicitn´ı metody Jednokrokov´e metody

y

3

[x3 , Y ]

[x3 , y(x3)]

y=y(x)

x x

I I

0

x

1

x

x

3

x

4

Pˇresn´ e ˇreˇsen´ı: y = y (x) pro x ∈ ha, bi. Pˇribliˇ zn´ e ˇreˇsen´ı Y i : krok h = (b − a)/n, xi = a + ih,

I

2

Y i ≈ y (xi )

Jednokrokov´ e metody: Y i+1 = Y i + hΦ(xi , Y i , h)

I I I I

Pˇr´ır˚ ustkov´a funkce Φ = Φ(x, y , h), pˇresn´y relativn´ı pˇr´ır˚ ustek ∆(x, y , h). Lok´ aln´ı diskretizaˇ cn´ı chyba: δ = Φ − ∆ Glob´ aln´ı (akumulovan´ a) chyba: e(xi ) = Y i − y (xi ) Konvergence, ˇr´ad metody, stabilita.

Pˇr. Eulerova metoda Pˇr´ıklad

22. Je d´ana Cauchyova u ´loha y 0 (y − 4) = x +

p 3 1 − y,

y (3) = 2,

a) S krokem h = 0.5 urˇcete pˇribliˇznou hodnotu y (2) pomoc´ı Eulerovy metody.

Pˇr. Eulerova metoda Pˇr´ıklad

22. Je d´ana u ´loha x¨ + 2x˙ + 26x = (1 − t)+ , x(0) = 0.001, x(0) = 0, a) Volte krok a vypoˇc´ıtejte pˇribliˇzn´e ˇreˇsen´ı v intervalu [0,4]. b) Zjistˇete maximum numerick´eho ˇreˇsen´ı a zkuste odhadnout jeho chybu (ˇreˇste na poˇc´ıtaˇci).

Eulerova metoda Pˇr´ıklad

20. Je d´ana Cauchyova u ´loha   y1 tgx + y3 − 2 ~y 0 =  y1 + y2 ln (x + 1)  1 y3 y1 + 2y2 − x−2



 −1 ~y (1) =  1  2

a) Ovˇeˇrte, ˇze Cauchyova u ´loha m´a pr´avˇe jedno ˇreˇsen´ı b) Zapiˇste interval I jej´ıho maxim´aln´ıho ˇreˇsen´ı c) S krokem h = 0.2 urˇcete pˇribliˇznou hodnotu ~y (1.2) pomoc´ı Eulerovy metody

Pˇredn´aˇska ˇc. 6

I

Princip jednokrokov´ych metod typu Runge-Kutty.

I

Collatzova metoda (E1), RK4. Konvergence metod.

I

Praktick´e uˇzit´ı metod.

Explicitn´ı a implicitn´ı Eulerova metoda. I

Explicitn´ı Eulerova metoda Y n+1 − Y n = f (xn , Y n ), h

I

Implicitn´ı Eulerova metoda Y n+1 − Y n = f (xn+1 , Y n+1 ), h

I

Y n+1 = Y n + hf (xn , Y n ),

Y n+1 −hf (xn+1 , Y n+1 ) = Y n

ˇ s´ıme u Stabilita: Reˇ ´lohu (λ ∈ C) y 0 = λy ,

y (0) = 1,

kde zn´ame analytick´e ˇreˇsen´ı. Jak se chov´a numerick´e ˇreˇsen´ı? I

Glob´ aln´ı chyba metody: O(h).

Eulerova expl. metoda a pˇresn´e ˇreˇsen´ı Rovnice x¨ + x = 0

I

Konvergence k ˇreˇsen´ı pro h = 2π/20, 2π/40, , 2π/80, 2π/160.

Eulerova expl. metoda, chyba Rovnice x¨ + x = 0

I

Konvergence k ˇreˇsen´ı pro h = 2π/20, 2π/40, , 2π/80, 2π/160.

Zobecnˇen´ı Eulerovy explicitn´ı metody Jednokrokov´e metody

y

3

[x3 , Y ]

[x3 , y(x3)]

y=y(x)

x x

I I

0

x

1

x

x

3

x

4

Pˇresn´ e ˇreˇsen´ı: y = y (x), x ∈ ha, bi. Pˇribliˇ zn´ e ˇreˇsen´ı Y i : xi = a + ih,

I

2

Y i ≈ y (xi )

Jednokrokov´ a metoda Y i+1 = Y i + hΦ(xi , Y i , h)

I I I I

Pˇr´ır˚ ustkov´a funkce Φ, pˇresn´y relativn´ı pˇr´ır˚ ustek ∆. Lok´ aln´ı diskretizaˇ cn´ı chyba δn , konzistentn´ı metoda. Glob´ aln´ı chyba (akumulovan´ a diskretizaˇ cn´ı chyba): ei = Y i − y (xi ) Konvergence, ˇr´ad metody, stabilita.

Zobecnˇen´ı Eulerovy explicitn´ı metody Metody vyˇsˇs´ıho ˇr´ adu

y

3

[x3 , Y ]

[x3 , y(x3)]

y=y(x)

x x

I I

0

x

1

x

2

x

3

x

4

Cauchyova u ´loha y 0 = f (x, y ), y (x0 ) = y0 Metody Taylorova rozvoje 1 y (x0 + h) = y (x0 ) + y 0 (x0 )h + y 00 (x0 )h2 + O(h3 ) 2

I

Zderivov´an´ım dostaneme d y 00 (x) = (f (x, y (x))) = fx (x, y (x)) + fy (x, y (x))y 0 (x) dx Tedy ˜ y , h) = f (x, y ) + hfx (x, y ) + hfy (x, y )f (x, y ) Φ(x,

Princip jednokrokov´ych metod typu Runge-Kutta: n = 2.

I

Metoda RK2 Y i+1 = Y i + h(ω1 k1 + ω2 k2 )

I

k1 = f (xi , Y i ) k2 = f (xi + α2 h, Y i + hβ21 k1 ) I

˜ 4 koeficienty: ω1 , ω2 , α2 , β21 . Srovn´ame Φ s Φ,

I

Jejich volbou dos´ahneme 2. ˇr´ad pˇresnosti (srovn´an´ım s metodou Taylorova rozvoje)

Princip jednokrokov´ych metod typu Runge-Kutty.

Collatzova metoda: Y i+1 = Y i + hk2 k2 = f (xi + h2 , Y i + h2 k1 ), k1 = f (xi , Y i ), I

grafick´e zn´azornˇen´ı,

I

ˇr´ad metody: O(h2 ).

Collatzova metoda - srovn´an´ı s pˇresn´ym ˇreˇsen´ım Rovnice x¨ + x = 0

I

Konvergence k ˇreˇsen´ı pro h = 2π/20, 2π/40, , 2π/80.

Collatzova metoda - chyby v z´avislosti na kroku h Rovnice x¨ + x = 0

I

Konvergence k ˇreˇsen´ı pro h = 2π/20, 2π/40, , 2π/80.

Princip jednokrokov´ych metod typu Runge-Kutta I

Jednokrokov´ a metoda n-t´ eho stupnˇ e: Y i+1 = Y i + hΦ(xi , Y i , h)

I

Metody RK: Speci´aln´ı volba pˇr´ır˚ ustkov´e funkce Φ(xi , Y i , h) i

Φ(xi , Y , h) =

n X

ωi ki

i=1 I

kde i

ki = f (xi + αi h, Y + h

i−1 X

βij ki )

j=1 I

Volba koeficient˚ u: co nejvyˇsˇs´ıho ˇr´ad konvergence, z´avis´ı na n.

Metoda Runge-Kutty 4. ˇr´adu.

Runge-Kutta 4. ˇr´ adu: h Y i+1 = Y i + (k1 + 2k2 + 2k3 + k4 ), 6 k1 = f (xi , Y i ), k2 = f (xi + h2 , Y i + h2 k1 ), k3 = f (xi + h2 , Y i + h2 k2 ), k4 = f (xi + h, Y i + hk3 ), ˇ ad metody: O(h4 ). R´

ˇ ad konvergence metod typu Runge-Kutty. R´ I

Y

i+1

i

=Y +h

n X

! ωi ki

,

i=1 I

kde ki = f (xi + αi h, Y i +

i−1 X

βij ki )

j=1

I

I

n p 1 1 2 2 ˇ ad konvergence: R´ 3 3 4 4 5 4 6 5 Pozn. Volba kroku: h lze volit i z´aporn´e, lze pouˇz´ıt i volba h > 1 - souvislost s O(hp )!

Pˇr´ıklady - Cauchyova u´loha. 21. Je d´ana Cauchyova u ´loha   2x√− y22 + y1 + 1 0 ~y = 4 − y1 − 2x

 ~y (0) =

0 1



a) Zapiˇste oblast G v n´ıˇz jsou splnˇeny postaˇcuj´ıc´ıpodm´ınky existence a jednoznaˇcnosti ˇreˇsen´ı Cauchyovy u ´lohy. b) Uˇzit´ım Collatzovy metody (1.modifikace Eulerovy metody) s krokem h = 1 urˇcete pˇribliˇznˇe hodnotu ˇreˇsen´ı v bodˇe x = 2. c)

Necht’ yi oznaˇ cuje hodnotu numerick´ eho ˇreˇsen´ı v bodˇ e xi z´ıskan´ eho Collatzovou metodou a yi∗ hodnotu pˇresn´ eho ˇreˇsen´ı. Zapiˇste pomoc´ı tˇ echto hodnot glob´ aln´ı chybu εi . Zapiˇste jak´ a je z´ avislost εi na kroku h a urˇ cete jak´ eho ˇr´ adu je Collatzova metoda. Odhadnˇ ete, jak se zmˇ en´ı glob´ aln´ı chyba v bodˇ e x pˇri zmˇ enˇ e kroku z h na h/4 u Collatzovy metody.

22. Je d´ana Cauchyova u ´loha y 000 − yy 00 +

x y0 ln x = 2 , y x −1

y (3) = −1, y 0 (3) = 0, y 00 (3) = 0.

Urˇcete s krokem h = 0, 4 pomoc´ı Collatzovy metody pˇribliˇznou hodnotu ˇreˇsen´ı y 00 (3.4).

Pˇredn´aˇska ˇc. 7

I

Metoda s´ıt´ı (koneˇcn´ych diferenc´ı) v 1 a 2D.

I

Okrajov´a u ´loha pro obyˇcejnou line´arn´ı diferenci´aln´ı rovnici 2. ˇr´adu. Existence a jednoznaˇcnost ˇreˇsen´ı.

I

Samoadjungovan´y tvar rovnice. Sturmovy okrajov´e podm´ınky.

Metoda s´ıt´ı 1D

a

b

x i

I

ˇ sen´ı hled´ame na I = ha, bi Reˇ Vol´ıme h = (b − a)/n, s´ıt’ xi = a + ih ˇ sen´ı aproximujeme y (xi ) ≈ Yi Reˇ

I

Derivace nahrad´ıme diferencemi, ˇreˇs´ıme soustavu

I I

~ ) = Fh Ah (Y

Metoda s´ıt´ı 2D

Ω

y x

I I

Aproximujeme u(x, y ), [x, y ] ∈ Ω ⊂ R2 : vol´ıme h > 0, oznaˇc´ıme xi = x0 + ih, s´ıt’ov´e ˇc´ary x = xi

Metoda s´ıt´ı 2D

Ω

y x

I

Aproximujeme u(x, y ), [x, y ] ∈ Ω ⊂ R2 : vol´ıme h > 0, oznaˇc´ıme xi = x0 + ih, s´ıt’ov´e ˇc´ary x = xi

I

oznaˇc´ıme yj = y0 + jh, s´ıt’ov´e ˇc´ary y = yj

I

Metoda s´ıt´ı 2D

Ω

y x

I I I I I

Aproximujeme u(x, y ), [x, y ] ∈ Ω ⊂ R2 : vol´ıme h > 0, oznaˇc´ıme xi = x0 + ih, s´ıt’ov´e ˇc´ary x = xi oznaˇc´ıme yj = y0 + jh, s´ıt’ov´e ˇc´ary y = yj oznaˇc´ıme Pi,j = [xi , yj ] pr˚ useˇc´ıky s´ıt’ov´ych ˇcar v Ω u(Pi,j ) ≈ Ui,j , derivace → diference, dostaneme ~ h (U~h ) = F ~h G

Okrajov´a u´loha pro diferenci´aln´ı rovnici 2. ˇr´adu Motivace: proˇc ˇreˇs´ıme?

I

Cauchyova u ´loha: my 00 = F (x, y , y 0 ),

y (a) = A, y 0 (a) = K

I

Fyzik´alnˇe: napˇr. hmotn´y bod [a, A], m´a danou rychlost, ?trajektorie?

I

Okrajov´ au ´loha: my 00 = F (x, y , y 0 ),

I

y (a) = A, y (b) = B

Fyzik´alnˇe: jak se z [a, A] dostat do [b, B] ?

Okrajov´a u´loha pro line´arn´ı diferenci´aln´ı rovnici 2. ˇr´adu Formulace I

Line´arn´ı diferenci´aln´ı rovnici 2. ˇr´adu y 00 + f1 (x)y 0 + f2 (x)y = g (x),

I

pro x ∈ (a, b),

Okrajov´e podm´ınky (Sturmovy okrajov´e podm´ınky), |α1 | + |α2 | = 6 0, |β1 | + |β2 | = 6 0,

α1 y 0 (a) + α2 y (a) = α3 , I

I

β1 y 0 (b) + β2 y (b) = β3 .

Speci´aln´ı volby: Okrajov´e podm´ınky Neumannovy

y 0 (a) = α,

y 0 (b) = β,

Dirichletovy

y (a) = α,

y (b) = β,

y (a) = α,

y 0 (b) = β.

Kdy existuje jednoznaˇcn´e ˇreˇsen´ı y (x)? pro okrajovou u ´lohu mnohem komplikovanˇ ejˇs´ı

Existence ˇreˇsen´ı pro okrajovou u´lohu Line´ arn´ı diferenci´ aln´ı rovnice 2. ˇr´ adu I

´ Uloha A: y 00 (x) = sin(x),

I

´ Uloha B: y 00 + y = 0,

I

y (0) = 0, y (π/2) = 1.

´ Uloha D: y 00 + y = 0,

I

y (0) = 0, y (π) = 1.

´ Uloha C: y 00 + y = 0,

I

y (0) = 2, y (π) = −2.

y (0) = 0, y (π) = 0.

´ Uloha E: −y 00 + y = 0,

y (0) = 0, y (π) = 0.

Ve vˇsech pˇr´ıpadech: hladk´a data. Existenci a jednoznaˇcnost ˇreˇsen´ı lze zaruˇcit jen pro u ´lohy ve speci´aln´ım tvaru.

Existence ˇreˇsen´ı pro okrajovou u´lohu Tvar rovnice v technick´ych u ´loh´ ach

V technick´ych u ´loh´ach m´a diferenci´aln´ı rovnice ˇcasto tvar −(py 0 )0 + qy = f ,

I

p, q - materi´alov´e konstanty, kladn´e/nez´aporn´e.

I

pro neizotropn´ı materi´al - p, q jako nez´aporn´e funkce p = p(x),

q = q(x).

Samoadjungovan´y tvar rovnice: −(p(x)y 0 )0 + q(x)y = f (x),

Existence a jednoznaˇcnost ˇreˇsen´ı pro u´lohy v samoadjungovan´em tvaru Je d´ana okrajov´a u ´loha: −(p(x)y 0 )0 + q(x)y = f (x), a okrajov´e podm´ınky α1 y 0 (a) + α2 y (a) = α3 ,

β1 y 0 (b) + β2 y (b) = β3 .

Necht’ plat´ı (i) p(x), p 0 (x), q(x), f (x) - spojit´e funkce na ha, bi (ii) p(x) > 0, q(x) ≥ 0 pro x ∈ ha, bi Pak existuje pr´avˇe jedno ˇreˇsen´ı (s v´yjimkou pˇr´ıpadu α2 = β2 = 0 ≡ q(x)).

Pˇrevod na samoadjungovan´y tvar Uvaˇ zujeme rovnici y 00 + f1 (x)y 0 + f2 (x)y = g (x),

pro x ∈ (a, b),

pˇren´asob´ıme −p(x) (nezn´am´a funkce) −p(x)y 00 −p(x)f1 (x)y 0 −p(x)f2 (x)y = −p(x)g (x), samoadjungovan´ y tvar −p(x)y 00 −p 0 (x)y 0 +q(x)y = f (x), {z } | −(py 0 )0

Dost´ av´ ame: −p 0 (x) = −p(x)f1 (x), q(x) = −p(x)f2 (x),

diferenci´aln´ı rovnice pro p(x) f (x) = −p(x)g (x).

Pˇrevod na samoadjungovan´y tvar. Pˇr´ıklady. 23. Formulujte Dirichletovu u ´lohu pro obyˇ cejnou line´ arn´ı diferenci´ aln´ı rovnici 2.ˇr´ adu v samoadjungovan´ em tvaru a) Uved’te podm´ınky pro jednoznaˇ cnou ˇreˇsitelnost u ´lohy (a) b) Zd˚ uvodnˇ ete, zda jsou postaˇ cuj´ıc´ı podm´ınky splnˇ eny pro u ´lohu 0 0

−(xy ) +

3−x x

y = − ln(2 + x)

y (1) = 0, y (2) = −4

24. Je d´ ana Dirichletova u ´loha y

00

−

2 0 2 y + (c − x)y + x = 0

x

y (2) = −1, y (4) = 3

a) Danou rovnici pˇreved’te na samoadjungovan´ y tvar b) Urˇ cete vˇsechny hodnoty parametru c ∈ R, pro nˇ eˇz jsou splnˇ eny postaˇ cuj´ıc´ı podm´ınky jednoznaˇ cn´ e ˇreˇsitelnosti Dirichletovy u ´lohy 25. Je d´ ana Dirichletova u ´loha y

00

−

4 x2 − 4

y (3) = 0, y (5) = −1

y =1

a) Ovˇ eˇrte, ˇze u ´loha m´ a pr´ avˇ e jedno ˇreˇsen´ı b) Odvod’te soustavu s´ıt’ov´ ych rovnic, kter´ a vznikne pˇri ˇreˇsen´ı dan´ eu ´lohy metodou s´ıt´ı s krokem h = 0.5 26. Je d´ ana rovnice y

00

+

2 0 y −

x

x 2+x

y =

1 x −3

a) Urˇ cete intervaly maxim´ aln´ıho ˇreˇsen´ı Cauchyovy u ´lohy pro danou rovnici b) Ovˇ eˇrte, ˇze jsou splnˇ eny postaˇ cuj´ıc´ı podm´ınky jednoznaˇ cn´ e ˇreˇsitelnosti Dirichletovy u ´lohy pro danou rovnici s okrajov´ ymi podm´ınkami y (−5) = −2, y (−3) = 0

Princip numerick´eho ˇreˇsen´ı −(p(x)y 0 )0 + q(x)y = f (x), a

+ okrajov´e podm´ınky b

x i

I I

Oznaˇcme pˇresn´ e ˇreˇsen´ı y = y (x) pro x ∈ I = ha, bi. Interval I rozdˇel´ıme ekvidistantnˇe s krokem h xi = a + i h,

I

oznaˇcme qi = q(xi ), fi = f (xi ), pi = p(xi )

Budeme hledat pˇribliˇ zn´ e ˇreˇsen´ı Y i ≈ y (xi )

I

Rovnici −(py 0 )0 + qy = f nahrad´ıme v bodech xi qy |xi ≈ qi Y i f |xi = fi − py


0 0

≈ dvoj´ım uˇzit´ım 1. centr´aln´ı diference

Diferenˇcn´ı n´ahrada v uzlu xi −(py 0 )0 + qy = f ,

x

x

i−1 I

x i−1/2

x i

x

i+1/2

i+1

1. centr´ aln´ı diference s krokem h/2 0 0

(py ) |xi ≈ y 0 (xi+ 1 ) ≈ 2

I

y (a) = A, y (b) = B,

pi+ 1 y 0 (xi+ 1 ) − pi− 1 y 0 (xi− 1 ) 2

2

2

h

Y i+1 − Y i , h

y 0 (xi− 1 ) ≈ 2

2

,

Y i − Y i−1 h

tedy dostaneme py


0 0

|xi ≈

pi− 1 Y i−1 − (pi+ 1 + pi− 1 )Y i + pi+ 1 Y i+1 2

2

2

h2

2

Princip numerick´eho ˇreˇsen´ı I

Oznaˇcme pˇresn´ e ˇreˇsen´ı y = y (x) pro x ∈ I = ha, bi.

I

Interval I rozdˇel´ıme ekvidistantnˇe s krokem h xi = a + i h,

I

oznaˇcme qi = q(xi ), fi = f (xi ), pi = p(xi )

Budeme hledat pˇribliˇ zn´ e ˇreˇsen´ı Y i ≈ y (xi )

I

Rovnici −(py 0 )0 + qy = f nahrad´ıme v bodech xi qy |xi ≈ qi Y i f |xi = fi −pi−1/2 Y i−1 + (pi+1/2 + pi−1/2 )Y i − pi+1/2 Y i+1 h2 Jak´a je lok´ aln´ı diskretizaˇ cn´ı chyba a glob´ aln´ı chyba? − py 0

I


0

≈

Pˇredn´aˇska ˇc. 8

I

Numerick´e ˇreˇsen´ı okrajov´e u ´lohy pro ODR.

I

Aproximace Dirichletov´ych a Neumannov´ych okrajov´ych podm´ınek, ˇ sen´ı vznikl´e soustavy line´arn´ıch rovnic. Reˇ

I

Princip numerick´eho ˇreˇsen´ı −(p(x)y 0 )0 + q(x)y = f (x), a

+ okrajov´e podm´ınky b

x i

I I

Oznaˇcme pˇresn´ e ˇreˇsen´ı y = y (x) pro x ∈ I = ha, bi. Interval I rozdˇel´ıme ekvidistantnˇe s krokem h xi = a + i h,

I

oznaˇcme qi = q(xi ), fi = f (xi ), pi = p(xi )

Budeme hledat pˇribliˇ zn´ e ˇreˇsen´ı Y i ≈ y (xi )

I

Rovnici −(py 0 )0 + qy = f nahrad´ıme v bodech xi qy |xi ≈ qi Y i f |xi = fi py


0 0

≈ dvoj´ım uˇzit´ım 1. centr´aln´ı diference s krokem h/2

Diferenˇcn´ı n´ahrada v uzlu xi −(py 0 )0 + qy = f ,

x

x

i−1 I

x i−1/2

x i

x

i+1/2

i+1

1. centr´ aln´ı diference s krokem h/2 0 0

(py ) |xi ≈ y 0 (xi+ 1 ) ≈ 2

I

y (a) = A, y (b) = B,

pi+ 1 y 0 (xi+ 1 ) − pi− 1 y 0 (xi− 1 ) 2

2

2

h

Y i+1 − Y i , h

y 0 (xi− 1 ) ≈ 2

2

,

Y i − Y i−1 h

tedy dostaneme py


0 0

|xi ≈

pi− 1 Y i−1 − (pi+ 1 + pi− 1 )Y i + pi+ 1 Y i+1 2

2

2

h2

2

Diferenˇcn´ı sch´ema pro Dirichletovu u´lohu Vlastnosti soustavy rovnic

Aproximovali jsme Dirichletovu u ´lohu

− (py 0 )0 + qy = f ,

y (a) = A, y (b) = B,

soustavou rovnic pro i = 1, . . . , n − 1

−pi− 1 Y i−1 + (pi+ 1 + pi− 1 + h2 qi )Y i − pi+ 1 Y i+1 = h2 fi 2

2

2

2

kde dosad´ıme za Y 0 = A, Y n = B. Vlastnosti soustavy rovnic: I

Matice soustavy je symetrick´a, tˇr´ıdiagon´aln´ı.

I

Matice soustavy je diagon´alnˇe dominantn´ı (p > 0, q ≥ 0).

I

Je-li q > 0 pak matice soustavy je ODD. Vˇzdy je IDD.

I

Matice soustavy je pozitivnˇe definitn´ı (p > 0, q ≥ 0).

(1)

N´ahrada Neumannovy okrajov´e podm´ınky Vlastnosti soustavy rovnic

Aproxinujeme u ´lohu

− (py 0 )0 + qy = f ,

I

y (a) = A, y 0 (b) = B,

(2)

diference 1. ˇr´ adu (ˇr´ad chyby!) Y n − Y n−1 =B h

I

1. centr´ aln´ı diference (ghost cell, pˇrid´ame xn+1 = b + h) Y n+1 − Y n−1 =B 2h −pn− 21 Y n−1 + (pn+ 12 + pn− 12 + h2 qn )Y n − pn+ 12 Y n+1 = h2 fn

I

I

3. bodov´a diference z xn , xn−1 , xn−2 (nevhodn´e: poruˇs´ıme DD i symetrii), v uzlu x = b postupujeme obdobnˇe.

Diferenˇcn´ı n´ahrada dif. rovnice

Vˇeta: Pro p ∈ C 3 (I ), y ∈ C 4 (I ) plat´ı 0 (py 0 )

= x=xi

pi+ 12 yi+1h−yi − pi− 12 yi −yh i−1 h

+ O(h2 )

Diferenˇcn´ı n´ahrada dif. rovnice Vˇeta: Pro p ∈ C 3 (I ), y ∈ C 4 (I ) plat´ı 0 (py 0 )

=

0

y (x

i+ 1 2

)=

h

x=xi

Dk. (py 0 )0

pi+ 12 yi+1h−yi − pi− 12 yi −yh i−1

=

pi+ 12 y 0 (xi+ 12 ) − pi− 12 y 0 (xi− 12 ) h

x=xi

yi+1 − yi h2 000 3 + O(h ) − y (xi+1/2 ), h 24

0

y (x

i− 1 2

)=

+ O(h2 )

+ O(h2 ),

yi − yi−1 h2 000 3 + O(h ) − y (xi−1/2 ), h 24

nebot’ plat´ı z(x0 + h/2) − z(x0 − h/2) = z 0 (x0 )h + 0

z(x0 + h/2) = z(x0 ) + z (x0 ) 0

z(x0 − h/2) = z(x0 ) − z (x0 )

h 2 h 2

2 h3 000 z (x0 ) + O(h4 ), 3! 8

+

1 00 h2 1 000 h3 4 z (x0 ) + z (x0 ) + O(h ) 2 4 3! 8

+

1 00 h2 1 000 h3 4 z (x0 ) − z (x0 ) + O(h ) 2 4 3! 8

Lok´aln´ı diskretizaˇcn´ı chyba, glob´aln´ı chyba Dirichletova u ´loha

V´ıme: Pˇresn´e ˇreˇsen´ı Dirichletovy u ´lohy y (x) splˇ nuje

− (py 0 )0 + qy = f ,

y (a) = A, y (b) = B,

pro y ∈ C 4 (I ) a ~y hodnoty ˇreˇsen´ı y (x) v uzlech s´ıtˇe xi plat´ı,

Ah~y = Fh + ~ηh ,

ηh,i = O(h2 )

~ splˇ Pˇribliˇzn´e ˇreˇsen´ı Y nuje ~ =F AY Lok´ aln´ı diskretizaˇ cn´ı chyba je O(h2 ). ~? Lze nˇ eco ˇr´ıct o glob´ aln´ı chybˇ e chybˇ e e = ~y − Y

(3)

Stabilita pro spec. rovnici I I

´ Uloha −y 00 = f (x), y (a) = 0, y (b) = 0. Vede na soustavu rovnic s matic´ı A 2 -1 0 0 0 0 0 0

I

-1 2 -1 0 0 0 0 0

0 0 -1 2 -1 0 0 0

0 0 0 -1 2 -1 0 0

0 0 0 0 -1 2 -1 0

0 0 0 0 0 -1 2 -1

0 0 0 0 0 0 -1 2

Uk´aˇzeme Au = λu, kde pro k = 1, . . . , n − 1 u = (uj )n−1 j=1 ,

I

0 -1 2 -1 0 0 0 0

uj = exp(±iπkj/n),

λ = 2−2 cos(πk/n)

Symetrick´a poz. definitn´ı matice: kA−1 k2 = 1/λMIN = 1/(2 − 2 cos(π/n)) ≈ n2 /π 2 = 1/(π 2 h2 )

Numerick´e ˇreˇsen´ı vybran´eho probl´emu

4 pi cos( x), 2 π 2 Zn´ame analytick´e ˇreˇsen´ı. −y 00 =

y (−1) = y (1) = 0

Numerick´a ˇreˇsen´ı, konvergence −y 00 =

4 pi cos( x), 2 π 2

y (−1) = y (1) = 0

Numerick´e ˇreˇsen´ı, krok h = 14 , h = 81 , h =

1 16 , h

=

1 32

Numerick´a ˇreˇsen´ı, chyba −y 00 =

4 pi cos( x), 2 π 2

y (−1) = y (1) = 0

Chyba numerick´eho ˇreˇsen´ı(log), krok h = 14 , h = 81 , h =

1 16 , h

=

1 32

Numerick´a ˇreˇsen´ı Konvergence z´ avis´ı na u ´loze

−(x 2 y 0 )0 = 1,

y (0) = y (1) = 0

Numerick´e ˇreˇsen´ı, krok h = 14 , h = 81 , h =

1 16 , h

=

1 32

Pˇr´ıklady. Metoda s´ıt´ı pro Dirichletovu u ´lohu v samoadjungovan´em tvaru 27. Je d´ ana Dirichletova u ´loha

y

00

−

2 0 2 y + (c − x)y + x = 0

x

y (2) = −1, y (4) = 3

a) Danou rovnici pˇreved’te na samoadjungovan´ y tvar b) Urˇ cete vˇsechny hodnoty parametru c ∈ R, pro nˇ eˇz jsou splnˇ eny postaˇ cuj´ıc´ı podm´ınky jednoznaˇ cn´ e ˇreˇsitelnosti Dirichletovy u ´lohy c) Napiˇste prvn´ı dvˇ e rovnice soustavy s´ıt’ov´ ych rovnic kter´ a vznikne pˇri ˇreˇsen´ı dan´ eu ´lohy metodou s´ıt´ı s krokem h = 0.2 pro c = 1 28. Je d´ ana Dirichletova u ´loha

y

00

−

4 x2 − 4

y (3) = 0, y (5) = −1

y =1

a) Ovˇ eˇrte, ˇze u ´loha m´ a pr´ avˇ e jedno ˇreˇsen´ı b) Odvod’te soustavu s´ıt’ov´ ych rovnic, kter´ a vznikne pˇri ˇreˇsen´ı dan´ eu ´lohy metodou s´ıt´ı s krokem h = 0.5 29. Je d´ ana rovnice y

00

+

2 0 y −

x

x 2+x

y =

1 x −3

a) Urˇ cete intervaly maxim´ aln´ıho ˇreˇsen´ı Cauchyovy u ´lohy pro danou rovnici b) Ovˇ eˇrte, ˇze jsou splnˇ eny postaˇ cuj´ıc´ı podm´ınky jednoznaˇ cn´ e ˇreˇsitelnosti Dirichletovy u ´lohy pro danou rovnici s okrajov´ ymi podm´ınkami y (−5) = −2, y (−3) = 0 c) Napiˇste prvn´ı dvˇ e rovnice soustavy s´ıt’ov´ ych rovnic kter´ a vznikne pˇri ˇreˇsen´ı dan´ eu ´lohy metodou s´ıt´ı s krokem h = 0.4

Pˇredn´aˇska ˇc. 9

I

Klasifikace line´arn´ıch parci´aln´ıch diferenci´aln´ıch rovnic 2.ˇr´adu dvou nez´avisle promˇenn´ych.

I

Numerick´e ˇreˇsen´ı Dirichletovy u ´lohy pro Poissonovu rovnici.

Klasifikace parci´aln´ıch diferenci´aln´ıch rovnic(PDR) Pˇr´ıklady PDR I

Eliptick´a rovnice: Poissonova rovnice (4u = f ) ∂2u ∂2u + =f, ∂x 2 ∂y 2 nebo −∇ · ∇u = f ,

I

Parabolick´a rovnice: Rovnice veden´ı tepla (p > 0) ∂u ∂2u =p 2 +f, ∂t ∂x

I

Hyperbolick´a rovnice: Vlnov´a rovnice (c > 0) 2 ∂2u 2∂ u = c +f, ∂2t ∂x 2

Dirichletova u´loha pro Poissonovu rovnici Existence ˇreˇsen´ı?

I

Hled´ame u ∈ C 2 (Ω) ∪ C 1 (Ω) takov´e, ˇze 4u = f v Ω,

u=ϕ

na ∂Ω

I

Existence (klasick´eho) ˇreˇsen´ı z´avis´ı na vlastnostech f , ϕ a Ω.

I

Vˇ eta: Ω oblast s hladkou hranic´ı (C 2 ), ϕ ∈ C (Ω), f ∈ C 1 (Ω), pak existuje pr´avˇe jedno ˇreˇsen´ı u ∈ C 2 (Ω).

I

Vˇ eta: Ω konvexn´ı, Lipschitzovsky spojit´a hranice, ϕ ∈ C (Ω), f ≡ 0, pak existuje pr´avˇe jedno ˇreˇsen´ı u ∈ C 2 (Ω).

Metoda s´ıt´ı 2D

Ω

y x

I I I I I

Aproximujeme u(x, y ), [x, y ] ∈ Ω ⊂ R2 : vol´ıme h > 0, oznaˇc´ıme xi = x0 + ih, s´ıt’ov´e ˇc´ary x = xi oznaˇc´ıme yj = y0 + jh, s´ıt’ov´e ˇc´ary y = yj oznaˇc´ıme Pi,j = [xi , yj ] pr˚ useˇc´ıky s´ıt’ov´ych ˇcar v Ω u(Pi,j ) ≈ Ui,j , derivace → diference, dostaneme ~ h (U~h ) = F ~h G

Diferenˇcn´ı n´ahrady parci´aln´ı derivace Pi,j+1 h h

Pi−1,j

Pi,j

Pi+1,j

Pi,j−1 I

2. centr´aln´ı diference f (x0 + h) − 2f (x) + f (x0 − h) = f 00 (x0 ) + O(h2 ). h2

I

Uˇzijeme 2. centr´aln´ı diferenci ve smˇeru x resp. y, chyba O(h2 ) nahrad´ıme

I

Ui−1,j − 2Ui,j + Ui+1,j ∂2u (xi , yj ) ≈ , 2 ∂x h2 Ui,j−1 − 2Ui,j + Ui,j+1 ∂2u (xi , yj ) ≈ , 2 ∂y h2

Diferenˇcn´ı n´ahrady 2. parci´aln´ı derivace Pi,j+1 h h

Pi−1,j

Pi,j

Pi+1,j

Pi,j−1 I

I

Ui−1,j − 2Ui,j + Ui+1,j ∂2u (xi , yj ) ≈ , ∂x 2 h2 Ui,j−1 − 2Ui,j + Ui,j+1 ∂2u (xi , yj ) ≈ , 2 ∂y h2 N´ahrada rovnice 4u = f v bodˇe Pi,j −Ui,j−1 − Ui−1,j + 4Ui,j − Ui,j+1 − Ui+1,j = −h2 fi,j

I

jen pro Pi,j , kter´y m´a vˇsechny sousedy“, takov´y uzel ” naz´yv´ame regul´ arn´ı.

Hraniˇcn´ı a neregul´arn´ı uzly

Ω

y x

I

Aproximujeme ˇreˇsen´ı u(Pi,j ) ≈ Ui,j

I

Oznaˇc´ıme hraniˇ cn´ı, regul´ arn´ı a neregul´ arn´ı uzly.

I

v hraniˇ cn´ıch uzlech Q - uˇzijeme Dirichletovy okrajov´e podm´ınky UQ = ϕ(Q)

I

v neregul´ arn´ıch uzlech PN - chceme zachovat ˇr´ad aproximace!

N´ahrada v neregul´arn´ım uzlu u

h R

I

δh N

Q

sitova primka

Hodnotu v neregul´arn´ım uzlu - line´arn´ı interpolace UR − UN UN − UQ = h δh

I

tedy (1 + δ)UN − δUR = UQ

I

Pro hladk´e ˇreˇsen´ı je tato n´ahrada 2. ˇr´adu pˇresnosti, tj. chyba O(h2 ). Dle Taylorova rozvoje v PN : u(x + δh) = u(x) + δhu 0 (x) + O(h2 ) u(x − h) = u(x) − hu 0 (x) + O(h2 )

ˇ ad aproximace Vlastnosti soustavy rovnic. R´ I

V´ysledn´a soustava rovnic je line´arn´ı soustava rovnic (n o n nezn´am´ych).

I

Matice soustavy je symetrick´a (?, pozor, neregul´arn´ı uzly)

I

Matice soustavy je DD (ne vˇsak ODD), vˇetˇsinou IDD.

I

Matice soustavy je pozitivnˇe definitn´ı.

I

Lok´aln´ı chyba aproximace O(h2 ). Co plat´ı pro glob´aln´ı chybu? eh = u − Uh Ah Uh = Fh Ah u = Fh + ηh ,

I

Je eh = O(h2 )?

ηh = O(h2 )

Poissonova rovnice Pˇr´ıklad

33. Je d´ana Dirichletova u ´loha ∆u = 4 v oblasti Ω tvoˇren´e ˇctyˇru ´heln´ıkem s vrcholy [0; 0], [2; 0], [0; 1.8], [1.4; 1.8] s okrajovou podm´ınkou u(x, y ) = x 2 na hranici Γ = ∂Ω b) Naˇcrtnˇete obr´azek s ˇc´ıslov´an´ım uzl˚ u pro h = 0.6 c) Sestavte s´ıt’ov´e rovnice v uzlech tak, aby metoda byla 2.ˇr´adu pˇresnosti 34. a) Je d´ana Dirichletova u ´loha ∆u = x(y + 1) v oblasti Ω tvoˇren´e ˇctyˇru ´heln´ıkem s vrcholy [0; 0], [1.8; 0], [0; 1.5] a [1.5; 1.5]. s okrajovou podm´ınkou u(x, y ) = x + y na hraniciΓ. (a) Volte h = 0.5, nakreslete obr´azek oblasti, zobrazte vˇsechny s´ıt’ov´e ˇc´ary, s´ıt’ov´e uzly uvnitˇr oblasti, regul´arn´ı neregul´arn´ı a hraniˇcn´ı uzly, ˇc´ıslov´an´ı uzl˚ u. (b) Sestavte s´ıt’ov´e rovnice, v neregul´arn´ıch uzlech uˇzijte line´arn´ı interpolace.

Pˇredn´aˇska ˇc. 10 Rovnice veden´ı tepla

I

Metoda s´ıt´ı pro rovnici veden´ı tepla.

I

Explicitn´ı a implicitn´ı sch´ema.

I

Konvergence a stabilita sch´emat.

Rovnice veden´ı tepla

I

´ Uloha:

I

poˇ c´ ateˇ cn´ı podm´ınka

∂u ∂2u =p 2 +f, ∂t ∂x

u(x, 0) = u0 (x), I

x ∈ ha, bi,

okrajov´ e podm´ınky u(a, t) = α(t),

u(b, t) = β(t),

t ≥ 0.

Metoda s´ıt´ı a aproximace ˇreˇsen´ı I

I

Rovnice:

∂2u ∂u =p 2 +f, ∂t ∂x S´ıt’: Vol´ıme krok h = (b − a)/n a krok τ > 0. xi = a + ih,

I

S´ıt’ov´e uzly Pik a aproximace ˇreˇsen´ı Uik Pik = [xi , tk ],

I

tk = kτ,

Uik ≈ u(xi , tk ).

Postup ˇreˇsen´ı: {U0i } → {Ui1 } → {Ui2 } → {Ui3 } . . .

Rovnice veden´ı tepla Explicitn´ı sch´ema

Pi k + 1 τ h k

Pi I

Rovnice:

I

N´ahrada v uzlu Pik , chyba O(h2 + τ )

∂2u ∂u =p 2 +f, ∂t ∂x

k − 2U k + U k Ui−1 ∂2u k i i+1 ) ≈ (P i ∂x 2 h2

U k+1 − Uik ∂u k (Pi ) ≈ i ∂t τ

Rovnice veden´ı tepla Explicitn´ı sch´ema

Pi k + 1 τ h k

Pi I

I

I

Rovnice:

∂u ∂2u =p 2 +f, ∂t ∂x k N´ahrada v uzlu Pi , chyba O(h2 + τ ) k U k − 2Uik + Ui+1 Uik+1 − Uik = p i−1 + fi k , τ h2 n´asob´ıme τ , oznaˇc´ıme σ = pτ h2 k k Uik+1 = σUi−1 + (1 − 2σ)Uik + σUi+1 + τ fi k ,

I

Stabilita: σ=

pτ 1 ≤ . 2 h 2

Rovnice veden´ı tepla Explicitn´ı sch´ema - stabilita I

Sch´ema: k k Uik+1 = σUi−1 + (1 − 2σ)Uik + σUi+1 + τ fi k ,

I

Je chyba ”mal´a”, je-li h a τ ”mal´e”? vol´ıme h = 0.01, τ = 0.001, σ = hτ2 = 10 ut = uxx ,

t x 0 0.001 0.002 0.003

0 0 0 0 0

2

u(x, 0) = (100x) ,

u(0, t) = 0,

0.01 1 21 -159 3461

0.03 9 29 49 ...

0.02 4 24 44 -1936

u(0.05, t) = 25

0.04 16 36 -144 ...4

0.05 25 25 25 25

Rovnice veden´ı tepla Explicitn´ı sch´ema - stabilita

I

Sch´ema: k k Uik+1 = σUi−1 + (1 − 2σ)Uik + σUi+1 + τ fi k ,

I

Je chyba ”mal´a”, je-li h a τ ”mal´e”? vol´ıme h = 0.01, τ = 0.00005, σ = hτ2 = 0.5 ut = uxx ,

t x 0 0.0005 0.0010 0.0015 0.0020

0 0 0 0 0 0

2

u(x, 0) = (100x) ,

0.01 1 2 2.5 3 3.375

0.02 4 5 6 6.75 7.625

u(0, t) = 0,

0.03 9 10 11 12.25 12.375

u(0.05, t) = 25

0.04 16 17 17.5 18 18.625

0.05 25 25 25 25 25

Rovnice veden´ı tepla Explicitn´ı sch´ema

Pi k + 1 τ h k

Pi I

Sch´ema: k k Uik+1 = σUi−1 + (1 − 2σ)Uik + σUi+1 + τ fi k ,

I

Chyba aproximace: O(h2 + τ ),

I

Stabilita σ=

1 pτ ≤ . 2 h 2

Rovnice veden´ı tepla Implicitn´ı sch´ema

h k+1

Pi τ

I

I

Rovnice:

∂u ∂2u =p 2 +f, ∂t ∂x k+1 N´ahrada v uzlu Pi , chyba O(h2 + τ ) k+1 k+1 Ui−1 − 2Uik+1 + Ui+1 ∂ 2 u k+1 (P ) ≈ ∂x 2 i h2

U k+1 − Uik ∂u k+1 (Pi ) ≈ i ∂t τ

Rovnice veden´ı tepla Implicitn´ı sch´ema

h k+1

Pi τ

I

Rovnice:

I

N´ahrada v uzlu Pik+1 , chyba O(h2 + τ )

∂2u ∂u =p 2 +f, ∂t ∂x

k+1 U k+1 − 2Uik+1 + Ui+1 Uik+1 − Uik = p i−1 + fi k+1 , τ h2 I

n´asob´ıme τ , oznaˇc´ıme σ =

pτ , h2

dost´av´ame

k+1 −σUi−1 + (1 + 2σ)Uik+1 − σUik+1 = Uik + τ fi k+1 ,

Rovnice veden´ı tepla Srovn´ an´ı explicitn´ı a implicitn´ı sch´ema

Pi k + 1

h k+1

Pi

τ h

τ

k

Pi I

Explicitn´ı sch´ ema: k k Uik+1 = σUi−1 + (1 − 2σ)Uik + σUi+1 + τ fi k ,

I I I

Chyba aproximace: O(h2 + τ ), ≤ 12 Stabilita: σ = pτ h2 Implicitn´ı sch´ ema: k+1 −σUi−1 + (1 + 2σ)Uik+1 − σUik+1 = Uik + τ fi k+1 ,

I I

Chyba aproximace: O(h2 + τ ), Stabilita: Vˇzdy.

Rovnice veden´ı tepla Srovn´ an´ı explicitn´ı a implicitn´ı sch´ema

Pi k + 1

h k+1

Pi

τ h

τ

k

Pi I

Explicitn´ı sch´ ema: k k Uik+1 = σUi−1 + (1 − 2σ)Uik + σUi+1 + τ fi k ,

I I I

Chyba aproximace: O(h2 + τ ), ≤ 12 Stabilita: σ = pτ h2 Implicitn´ı sch´ ema: k+1 −σUi−1 + (1 + 2σ)Uik+1 − σUik+1 = Uik + τ fi k+1 ,

I I I

Chyba aproximace: O(h2 + τ ), Stabilita: Vˇzdy. Lze Crank-Nicholson: Chyba aproximace: O(h2 + τ 2 ), stabilita vˇzdy.

Rovnice veden´ı tepla Pˇr´ıklad 35.a) Je d´ ana sm´ıˇsen´ au ´loha ∂u ∂t

= 0.3

∂2u ∂x 2

+ x + 2t

u(x, 0) = x

2

v oblasti Ω = {[x; t] : x ∈ (0; 1); t > 0}

pro x ∈ h0; 1i

u(0, t) = arctg(t),

u(1, t) =

1 2t + 1

pro t ≥ 0

b) Ovˇ eˇrte splnˇ en´ı podm´ınek souhlasu c) Urˇ cete τ a minim´ aln´ı krok h tak, aby pˇri jej´ım ˇreˇsen´ı stabiln´ı explicitn´ı metodou leˇzel bod P = [0.25; 0.1] v prv´ eˇ casov´ e vrstvˇ e d) Pro hodnoty τ a h z bodu (b) urˇ cete pˇribliˇznou hodnotu ˇreˇsen´ı v bodˇ e P uˇzit´ım explicitn´ı metody e) Pˇri h = τ = 0.25 sestavte soustavu s´ıt’ov´ ych rovnic pro prvn´ı ˇ casovou vrstvu uˇzit´ım implicitn´ı formule 36. Je d´ ana rovnice ∂u ∂t

= 2.5

∂2u ∂x 2

v oblasti Ω = {[x; t] : x ∈ (0; 2); t > 0}

b) Pˇri zadan´ ych podm´ınk´ ach u(x, 0) = x(2 − x),

u(0, t) = 30t,

u(2, t) = 0,

pro x ∈ h0; 2i, t ≥ 0,

sestavte soustavu s´ıt’ov´ ych rovnic pro prvn´ı ˇ casovou vrstvu pomoc´ı implicitn´ıho schematu. Volte h = 0.5 a τ = 0.1 c) Rozhodnˇ ete, zda lze volit ˇ casov´ y krok τ = 0.01, resp. τ = 1 aby pro dan´ y krok v ose x bylo uˇzit´ e schema stabiln´ı

Pˇredn´aˇska ˇc. 11

I

Metoda s´ıt´ı pro vlnovou rovnici.

I

Explicitn´ı a implicitn´ı sch´ema.

I

Konvergence a stabilita sch´emat.

Formulace sm´ıˇsen´e u´lohy pro vlnovou rovnici I

Vlnov´a rovnice (c > 0) 2 ∂2u 2∂ u = c +f, ∂2t ∂x 2

I

poˇ c´ ateˇ cn´ı podm´ınky u(x, 0) = ϕ(x),

I

∂u (x, 0) = ψ(x), ∂t okrajov´ e podm´ınky u(a, t) = α(t),

I

x ∈ ha, bi, x ∈ ha, bi,

u(b, t) = β(t),

Jsou splnˇeny podm´ınky souhlasu.

t ≥ 0.

Explicitn´ı schema P

k+1

i

P

k

k

P

i−1

P

i

k i+1

k−1

P i

I

Aproximace derivac´ı

I

2 ∂2u 2∂ u = c +f, ∂2t ∂x 2 S´ıt’: Vol´ıme krok h = (b − a)/n a krok τ > 0.

xi = a + ih,

tk = kτ,

Explicitn´ı sch´ema N´ ahrada v regul´ arn´ım uzlu

P

k+1

i

P

k

k

P

i−1

P

i

k i+1

k−1

P i

I

Rovnice:

I

N´ahrada v uzlu Pik :

2 ∂2u 2∂ u = c +f, ∂2t ∂x 2

Uk+1 − 2Uik + Uik−1 i

= c2

k − 2U k + U k Ui+1 i i−1

+fk

Explicitn´ı sch´ema N´ ahrada v regul´ arn´ım uzlu

P

k+1

i

P

k

i−1

k

P

P

i

k i+1

k−1

P i

I

N´ahrada v Uik+1 =

I

uzlu Pik , explicitn´ı sch´ ema 2 k 2 k k σ Ui+1 + 2(1 − σ )Ui + σ 2 Ui−1

Stabilita: σ=

cτ ≤ 1. h

− Uik−1 + τ 2 fi k

N´ahrada na 1. ˇcasov´e vrstvˇe N´ ahrada s chybou O(τ 2 )

I

N´ahrada s chybou O(τ ): U 1 − Ui0 u(xi , τ ) − u(xi , 0) ∂u (xi , 0) = + O(τ ) ≈ i ∂t τ τ

I

nebo-li Ui1 ≈ u(xi , τ ) = u(xi , 0) + τ

∂u (xi , 0) + O(τ 2 ) ∂t

Vlnov´a rovnice, explicitn´ı sch´ema

I

Lok´aln´ı chyba aproximace - 2. centr´aln´ı diference (viz Pˇredn´aˇska 7) O(h2 + τ 2 )

I

D˚ uleˇzit´a stabilita: σ ≤ 1

I

Lze uˇ z´ıt implicitn´ıho sch´ ematu, tj. aby byla zaruˇcena stabilita pro libovoln´e σ ≥ 0?

Implicitn´ı sch´ema N´ ahrada v regul´ arn´ım uzlu I

Rovnice:

I

N´ahrada v uzlu Pik :

2 ∂2u 2∂ u = c +f, ∂2t ∂x 2

Uik+1 − 2Uik + Uik−1 ∂2u k (P ) ≈ ∂2t i τ2 k−1 k−1 k−1 k+1 k+1 k+1 + Ui−1 + Ui−1 ∂2u k 1 Ui+1 − 2Ui 1 Ui+1 − 2Ui (P ) ≈ + ∂x 2 i 2 h2 2 h2

I

Dosad´ıme do rovnice, n´asob´ıme τ 2 , oznaˇc´ıme c 2τ 2 σ2 = 2 h

I

vyj´adˇr´ıme soustavu pro Uik+1

I

1 2 k+1 1 2 k+1 1 2 k−1 1 2 k−1 2 k+1 2 k−1 k 2 k − σ Ui−1 + (1 + σ )Ui − σ Ui+1 = σ Ui−1 − (1 + σ )Ui + σ Ui+1 + 2Ui + τ fi 2 2 2 2

Vlnov´a rovnice Pˇr´ıklad 37. Je d´ ana sm´ıˇsen´ au ´loha

∂2u ∂t 2

2

u(x, 0)

=

x ,

u(−1, t)

=

1,

=4

∂2u ∂x 2

+ x sin t

∂u

2

(x, 0) = 1 − x pro x ∈ h−1; 1i ∂t u(1, t) = cos t pro t ∈ h0; ∞)

Ovˇ eˇrte splnˇ en´ı podm´ınek souhlasu (pro polohu a rychlost) b) Urˇ cete maxim´ aln´ı krok τ tak, aby byla splnˇ ena podm´ınka stability pro explicitn´ı metodu s prostorov´ ym krokem h = 0.2 c) Odvod’te s´ıt’ov´ e rovnice pro prvn´ı ˇ casovou vrstvu pˇri n´ ahradˇ e ∂u (x, 0) s chybou O(τ ) ∂t d) Stanovte pˇribliˇznou hodnotu ˇreˇsen´ı v bodˇ e A = [0.2; 0.2]. Uˇzijte v´ ysledky z bod˚ u (b) a (c). 38.a) Je d´ ana sm´ıˇsen´ au ´loha ∂2u ∂t 2

u(x, 0)

=

x(x − 1) ,

u(0, t)

=

sin t ,

=4

∂u

∂2u ∂x 2

+x

2

(x, 0) = (1 − x) pro x ∈ h0; 1i ∂t u(1, t) = 0 pro t ∈ h0; ∞)

Ovˇ eˇrte splnˇ en´ı podm´ınek souhlasu. b) Pro explicitn´ı metodu volte h = 0.2. Urˇ cete τ tak, aby byla splnˇ ena podm´ınka stability a bod A = [0.4; 0.2] byl uzlem s´ıtˇ e c) Stanovte pˇribliˇznou hodnotu ˇreˇsen´ı v bodˇ e A. Pro prvn´ı ˇ casovou vrstvu uˇzijte n´ ahradu s chybou O(τ ).

Pˇredn´aˇska ˇc. 12

I

Aproximace metodou nejmenˇs´ıch ˇctverc˚ u

Aproximace pomoc´ı metody nejmenˇs´ıch ˇctverc˚ u 1

0.5

0

-0.5

-1

-1.5

-2

-2.5 0

I I

0.5

1

1.5

2

´ Uloha: Chceme naj´ıt z´avislost v namˇeˇren´ych datech obsahuj´ı nepˇresnosti, Princip aproximace: Hled´ame takovou funkci v dan´e(m) mnoˇzinˇe(prostoru), kter´a minimalizuje odchylku.

Aproximace pomoc´ı metody nejmenˇs´ıch ˇctverc˚ u 1

0.5

0

-0.5

-1

-1.5

-2

-2.5 0

0.5

1

1.5

I

Princip: Minimalizace kvadratick´e odchyklky X (p(xi ) − yi )2 δ 2 (p(x)) =

I

Volba funkce: Napˇr. polynom stupnˇe n

i

2

Aproximace pomoc´ı metody nejmenˇs´ıch ˇctverc˚ u I

Princip: Minimalizace kvadratick´e odchyklky X δ 2 (p(x)) = (p(xi ) − yi )2 i

I

d´ana tabulka dat [xi , yi ], minimalizujeme kvadratickou odchylku X G (a0 , a1 ) := δ 2 (p(x)) = (p(xi ) − yi )2 i

I I

odvozen´ı norm´aln´ıch rovnic ∂G /∂ak = 0 pro k = 0, 1. soustava norm´aln´ıch rovnic X X X a0 ( 1) + a1 ( xi ) = yi , i

a0 (

X i

i

i

X X xi ) + a1 ( xi2 ) = xi yi , i

i

Aproximace pomoc´ı metody nejmenˇs´ıch ˇctverc˚ u 13.a) Vysvˇetlete princip metody nejmenˇs´ıch ˇctverc˚ u pˇri aproximaci tabulky hodnot (xi , yi ) pro i = 1, . . . , n polynomem nejv´yˇse 1. stupnˇe. b) Odvod’te obecnˇe soustavu norm´aln´ıch rovnic pro pˇr´ıpad (a). c) Urˇcete polynom p1∗ nejv´yˇse 1. stupnˇe, kter´y ve smyslu metody nejmenˇs´ıch ˇctverc˚ u nejl´epe aproximuje danou tabulku hodnot: xi yi

-1 0.5

-1 -0.4

0 0.7

1 0.5

1 0.5

2 -0.4

14.a) Vysvˇetlete princip metody nejmenˇs´ıch ˇctverc˚ u pˇri aproximaci tabulky hodnot (xi , yi ) pro i = 1, . . . , n polynomem nejv´yˇse 2. stupnˇe. b) Odvod’te obecnˇe soustavu norm´aln´ıch rovnic pro pˇr´ıpad (a). c) Urˇcete polynom p1∗ nejv´yˇse 2. stupnˇe, kter´y ve smyslu metody nejmenˇs´ıch ˇctverc˚ u nejl´epe aproximuje danou tabulku hodnot. Urˇcete odpov´ıdaj´ıc´ı kvadratickou odchylku. xi yi

-2 9.9

-1 4

-1 4.1

0 0.1

0 0.2

1 -2

1 -2.5

2 -1.8

Opakov´an´ı

I

Pˇr´ıklady.

P´ısemka (1., 2.) 1. Je d´ ana tabulka hodnot xi yi

−2 0.8

−1 −0.58

−1 −0.62

0 −1.2

0 −0.5

0 −1.3

1 0

1 −0.8

2 1.2

a) Zapiˇste podm´ınku, kterou m´ a splˇ novat polynom nejv´ yˇse 2. stupnˇ e, kter´ y aproximuje tabulku hodnot xi yi metodou nejmenˇs´ıch ˇ ctverc˚ u. Odvod’te soustavu norm´ aln´ıch rovnic pro tento pˇr´ıpad. [8b] b) Sestavte soustavu norm´ aln´ıch rovnic pro zadanou tabulku hodnot a pro aproximaci polynomem p2∗ (x) nejv´ yˇse 2. stupnˇ e. Proved’te LU rozklad matice soustavy a uˇzit´ım LU rozkladu vypoˇ ctˇ ete pˇresn´ e ˇreˇsen´ı soustavy. [10b] c) Urˇ cete polynom p2∗ (x) nejv´ yˇse 2. stupnˇ e, kter´ y danou tabulku hodnot aproximuje nejl´ epe ve smyslu metody nejmenˇs´ıch ˇ ctverc˚ u. [7b] 2. Je d´ ana Cauchyova u ´loha 0 ~ y =



2x√− y22 + y1 + 1 4 − y1 − 2x



 ~ y (0) =

0 1



a) Zapiˇste oblast G v n´ıˇz jsou splnˇ eny postaˇ cuj´ıc´ıpodm´ınky existence a jednoznaˇ cnosti ˇreˇsen´ı Cauchyovy u ´lohy [5b] b) Uˇzit´ım Collatzovy metody (1.modifikace Eulerovy metody) s krokem h = 1 urˇ cete pˇribliˇznˇ e hodnotu ˇreˇsen´ı v bodˇ e x = 2. [10b] c) Necht’ yi oznaˇ cuje hodnotu numerick´ eho ˇreˇsen´ı v bodˇ e xi z´ıskan´ eho Collatzovou metodou a yi∗ hodnotu pˇresn´ eho ˇreˇsen´ı. Zapiˇste pomoc´ı tˇ echto hodnot glob´ aln´ı chybu εi . Zapiˇste jak´ a je z´ avislost εi na kroku h a urˇ cete jak´ eho ˇr´ adu je Collatzova metoda. Odhadnˇ ete, jak se zmˇ en´ı glob´ aln´ı chyba v bodˇ e x pˇri zmˇ enˇ e kroku z h na h/4 u Collatzovy metody. [10b]

P´ısemka (3., 4.) 3. Je d´ ana Dirichletova okrajov´ au ´loha

y

00

+

1 2x

2 0 y − 4y = − , x

y (1) = 1, y (5) = 0

a) Danou rovnici pˇrevedte na samoadjungovan´ y tvar. Zapiˇste postaˇ cuj´ıc´ı podm´ınky pro existenci a jednoznaˇ cnost ˇreˇsen´ı okrajov´ eu ´lohy. Ovˇ eˇrte zda jsou splnˇ eny. [8b] y (x+h)−y (x−h)

b) Ukaˇzte, ˇze v´ yraz je aproximac´ı y 0 (x) pro dostateˇ cnˇ e hladkou funkci y (x). Uˇzit´ım 2h tohoto vztahu odvod’te tvar s´ıt’ov´ ych rovnic pro diskretizaci u ´lohy v samoadjungovan´ em tvaru metodou s´ıt´ı s krokem h. [8b] c) Volte h = 1 a sestavte s´ıt’ov´ e rovnice pro Dirichletovu u ´lohu v samoadjungovan´ em tvaru z´ıskanou v a). Je pro tuto soustavu rovnic Gauss-Seidelova iteraˇ cn´ı metoda konvergentn´ı? Zd˚ uvodnˇ ete. [9b] 4. D´ ana sm´ıˇsen´ au ´loha ∂u ∂t u(x, 0) = ax + b u(0, t) = 2 + t

=2

∂2u ∂x 2

+

2x t+1

u(4, t) = −2

,

pro x ∈< 0, 4 >, pro t ∈ h0, 10i

a) Urˇ cete, pro kter´ e hodnoty a, b ∈ R jsou splnˇ eny podm´ınky souhlasu a odvod’te s´ıt’ovou rovnici v regul´ arn´ım uzlu (k + 1)-n´ı ˇ casov´ e vrstvy pˇri ˇreˇsen´ı dan´ eu ´lohy implicitn´ı metodou. Zapiˇste podm´ınku stability pro implicitn´ı metodu. [10b] b) Volte krok h a τ maxim´ aln´ı tak, aby bod A = [1; 21 ] byl uzlem s´ıtˇ e implicitn´ı metoda byla stabiln´ı. Sestavte rovnice implicitn´ıho schematu na prvn´ı ˇ casov´ e vrstvˇ e. [9b] c) Pro z´ıskanou soustavu rovnic z b) proved’te jeden krok Jacobiho iteraˇ cn´ı metody, poˇ c´ ateˇ cn´ı iteraci volte jako nulov´ y vektor (tj. X (0) = ~ 0). [6b]