ˇ MA4 Pˇredmet: Dnešní látka: ◮
Nehomogenní okrajové podmínky.
◮
Pokraˇcování OÚ pro PDR (jen pro fajnšmekry).
◮
ˇ ˇrešení. Jednoznaˇcnost zobecneného Metoda sítí v 1D.
◮
◮ ◮
Pˇribližné ˇrešení okrajových úloh. Aproximace vlastních cˇ ísel.
Literatura: ◮
Kapitoly 3, 4 a 2 d) ze skript Karel Rektorys: Matematika ˇ 43, CVUT, Praha, 2001.
◮
Text pˇrednášky na webové stránce pˇrednášejícího.
Obecná Dirichletova okrajová podmínka (2D, 3D) Klasické ˇrešení úlohy: najít u ∈ C 2 (Ω) ∩ C(Ω), aby −∆u = f v Ω ⊂ R2 , u = g na Γ.
ˇ Slabé ˇrešení úlohy: Vezmeme takovou funkci w ∈ W21 (Ω), aby ˇ splnovala w = g na Γ, pak hledáme u ∈ W21 (Ω), pro niž Z
Ω
˚ 1 (Ω), u−w ∈W Z2 ∇u · ∇v dx = fv dx Ω
˚ 1 (Ω). ∀v ∈ W 2
˚ 1 (Ω): Lze ˇrešit napˇr. užitím u = w + u0 , kde u0 ∈ W 2 Z Z Z ˚ 1 (Ω). ∇u0 · ∇v dx = fv dx − ∇w · ∇v dx ∀v ∈ W 2 Ω
Ω
Ω
Pˇríklad v 1D:
Okrajová úloha −u ′′ + ex u = cos x,
u(0) = 1, u(3) = −5.
Zvolíme napˇríklad w (x) = 1 − 2x, tj. u = u0 + w , navíc dokonce u ′′ = u0′′ . Pak OÚ pro neznámou funkci u0 : −u0′′ + ex (u0 + w ) = cos x,
u0 (0) = 0, u0 (3) = 0.
˚ 1 ([0, 3]), aby ∀v ∈ W ˚ 1 ([0, 3]) Slabá formulace: Najít u0 ∈ W 2 2 Z
3 0
u0′ (x)v ′ (x) +
x
e u0 (x)v(x) dx =
Z
3
(cos x−ex (1−2x))v(x) dx.
0
ˇ Rešení puvodní ˚ úlohy je u = u0 + 1 − 2x ∈ W21 ([0, 3]).
Inspirace i pro 1D OÚ s jinými OP: ′ 3 ′ − (2 + x )u + u = cos πx,
u(0) = 7, u ′ (1) = 4. Zvolme w (x) = 7 + 4x, tj. u = u0 + w a platí (2). Pak z ′ − (2 + x 3 )(u0 + w )′ + u0 + w = cos πx, ′ ′ − (2 + x 3 )u0′ − (2 + x 3 )w ′ + u0 + w = cos πx, ′ − (2 + x 3 )u0′ − 12x 2 + u0 + 7 + 4x = cos πx,
odvodíme OÚ pro neznámou funkci u0 : ′ − (2 + x 3 )u0′ + u0 = cos πx + 12x 2 − 4x − 7, u0 (0) = 0, u0′ (1) = 0.
Pˇríslušný operátor je sym. a poz. def., OÚ (1)-(2) ˇrešíme ˇ ˇ Rešení standardne. puvodní ˚ úlohy je u = u0 + 7 + 4x.
(1) (2)
Jiné okrajové podmínky ve 2D OÚ (pro fajnšpekry) Neumannova okrajová podmínka ∂u = h na Γ. ∂ν Slabou formulaci úlohy −∆u = f v Ω, ∂u ∂ν = h na Γ lze odvodit postupem naznaˇceným v minulé pˇrednášce (Jiná ˇ ˇrešení) a uplatnením ˇ ˇ cesta k zobecnenému Greenovy vety (viz též oddíl Newtonova okrajová podmínka na další stránce). ˇ Dospejeme k úloze: Najít funkci u ∈ HA = W21 (Ω) takovou, aby platilo Z Z Z ∂u ∂v ∂u ∂v + dx = fv dx + hv dS ∂x2 ∂x2 Ω ∂x1 ∂x1 Ω Γ ˇ Pro podmínka R existenci R ˇrešení je nutná doplnující f dx + h dS = 0. Ω Γ
∀v ∈ HA .
Newtonova (Robinova) okrajová podmínka ∂u + αu = q na Γ. ∂ν Slabou formulaci úlohy −∆u = f v Ω, ∂u ∂ν + αu = q na Γ lze odvodit postupem naznaˇceným v minulé pˇrednášce (Jiná ˇ ˇrešení) a uplatnením ˇ ˇ cesta k zobecnenému Greenovy vety. Vynásobme tedy obeˇ strany rovnice −∆u = f testovací funkcí v a integrujme R pˇres oblast Ω. Levou stranu, tj. − Ω ∆uv dx upravíme pomocí Greenovy ˇ vety:
−
Z 2 ∂ u ∂ 2u ∆uv dx = − v + v dx ∂x12 ∂x22 Ω Ω Z Z ∂u ∂u ∂u ∂v ∂u ∂v =− v ν1 + v ν2 dS + + dx ∂x1 ∂x2 ∂x1 ∂x1 ∂x2 ∂x2 Γ Ω Z Z ∂u =− v dS + ∇u · ∇v dx ∂ν Ω ZΓ Z = − (q − αu) v dS + ∇u · ∇v dx (viz okraj. podm.).
Z
Γ
Ω
R ˇ Clen Γ qv dS pˇrevedeme na druhou stranu rovnice, tj. k Ω fv dx . Slabá formulace: Najít funkci u ∈ W21 (Ω) takovou, aby platilo Z Z Z Z ∇u · ∇v dx + αuv dS = fv dx + qv dS ∀v ∈ W21 (Ω). R
Ω
Γ
Ω
Γ
Pokud α > c > 0, je pˇríslušný operátor pozitivneˇ definitní.
Z
Ω
∇u · ∇v dx +
Z
αuv dS = Γ
Z
Ω
fv dx +
Z
Γ
qv dS
∀v ∈ W21 (Ω).
Tato okrajová úloha modeluje napˇríklad ustálené vedení tepla v homogenním a izotropním materiálu s jednotkovou tepelnou vodivostí a s koeficientem α pˇrestupu tepla do ˇ ˇ vnejšího prostˇredí, u je teplota v telese Ω.
Kdy existuje slabé ˇrešení? Operátor A pˇríslušný této úloze je symetrický a pozitivneˇ definitní. Z Z (u, v )A = ∇u · ∇v dx + αuv dS Ω
Γ
Symetrie zˇrejmá. Pozitivní definitnost (s využitím Friedrichsovy nerovnosti), dokonce vzhledem k normeˇ k · kW21 (Ω) : 1 (u, u)A ≥ min( , α) 2
Z
Z
2
1 + 2
∇u · ∇u dx + u dS Z Γ 1 2 ≥ c min(1/2, α)kukL2 (Ω) + ∇u · ∇u dx 2 Ω ≥ cˆkuk2W 1(Ω) , Ω
Z
Ω
∇u · ∇u dx
2
kde c je konstanta z Friedr. nerovnosti (viz minulou pˇrednášku) a cˆ = min(c/2, cα, 1/2) > 0 je konstanta.
Proˇc se používá slabá formulace? ◮ Obecnejší ˇ úlohy než pˇri minimalizaci funkcionálu energie (pro symetrické pozitivneˇ definitní operátory jsou však možné obeˇ cesty a vedou k témuž ˇrešení). ◮
Existence ˇrešení za dosti obecných podmínek (Laxovo-Milgramovo lemma). Jednoznaˇcnost ˇrešení.
◮
ˇ eˇ pruhledné" "Pomern ˚ odvození z klasické formulace.
◮
Velmi vhodná pro teoretickou analýzu (odhady chyby, rychlost konvergence).
◮
Teoretický základ metody koneˇcných prvku. ˚
ˇ ˇrešení Jednoznaˇcnost slabého (zobecneného) ˇ Necht’ operátor A je pozitivneˇ definitní na prostoru V . Mejme dveˇ ˇrešení u1 , u2 ∈ V , tj. (u1 , v)A = (f , v) (u2 , v)A = (f , v)
∀v ∈ V , ∀v ∈ V .
Po odeˇctení (u1 − u2 , v)A = 0 ∀v ∈ V . ˇ Vezmeme v = u1 − u2 (jest v ∈ V ), pak 0 = (u1 − u2 , u1 − u2 )A ≥ cku1 − u2 k2A =⇒ u1 = u2 .
ˇ úlohy se smíšenými okrajovými podmínkami Ukázka složitejší 2 X ∂ ∂u aij + bu = f v Ω ⊂ R2 , − ∂xi ∂xj i,j=1
u = g na Γ1 ,
2 X
aij
i,j=1
∂u νi = h na Γ2 , ∂xj
2 X
i,j=1
aij
∂u νi + αu = q na Γ3 , ∂xj
kde Γ = Γ1 ∪ Γ2 ∪ Γ3 a meas1 Γi > 0, i = 1, 2, 3, dále aij > 0, b ≥ 0. Slabá formulace: Najdi funkci u ∈ W21 (Ω): u − w ∈ V , kde w ∈ W21 (Ω), w = g na Γ1 ; Z X Z Z 2 ∂u ∂v aij dx + buv dx + αuv dS ∂xj ∂xi Ω i,j=1 Ω Γ3 Z Z Z = fu dx + hv dS + qv dS ∀v ∈ V , Ω
Γ2
Γ3
˚ 1 (Ω) ⊂ V ⊂ W 1 (Ω). kde V = v ∈ W21 (Ω) : v |Γ1 = 0 . Jest W 2 2
Metoda sítí pro pˇribližné rˇ ešení 1D OÚ a(x)y ′′ (x)+b(x)y ′ (x)−q(x)y(x) = f (x),
x ∈ [a, b], a OP (3)
Pˇredpoklady: a, b, q, f ∈ C([a, b]), y ∈ C 4 ([a, b]) (kvuli ˚ odhadum ˚ chyby), a(x ) ≥ a0 > 0, q(x ) ≥ 0 ∀x ∈ [a, b].
ˇ verzí rovnice Poznámka: Rovnice (3) je (až na znaménko) obecnejší ′ ′ −(p(x )y (x )) + q(x )y (x ) = f (x ), p(x ) ≥ p0 > 0, q(x ) ≥ 0.
Myšlenka: ◮
◮
ˇ interval [a, b], Ekvidistantneˇ rozdelit a = x0 < x1 < · · · < xN = b, xi = a + ih, h = (b − a)/N. ˇ rovnice (1) jen v bodech xi (neužívá se Požadovat splnení tedy variaˇcní formulace!), a to jen pˇribližneˇ .
◮
Místo y ′ (xi ), y ′′ (xi ) použít pˇribližnou derivaci vyjádˇrenou diferenˇcním podílem.
◮
Sestavit lineární algebraické rovnice pro neznámé Yi , které aproximují yi ≡ y(xi ).
◮
Vyˇrešit výslednou soustavu lineárních algeb. rovnic.
Tayloruv ˚ polynom 2.5 2 1.5 f pol. 1 pol. 2 pol. 3 pol. 4 pol. 5
1 0.5 0 −0.5 −1 −1.5 −2
−0.5
0
0.5
1
1.5
Rozvoj funkce f v okolí bodu a: n X f (k ) (a) f (x) ≈ (x − a)k k! k =0
2
Aproximace derivací pomocí Taylorova rozvoje yi+1 = yi + hy ′ (xi ) +
h2 ′′ h3 h4 y (xi ) + y ′′′ (xi ) + y (4) (xi + θi+ h), 2! 3! 4! (4)
yi−1 = yi − hy ′ (xi ) +
h2 ′′ h3 h4 y (xi ) − y ′′′ (xi ) + y (4) (xi − θi− h), 2! 3! 4! (5)
kde θi+ ∈ (0, 1), θi− ∈ (0, 1). − yi−1 y (4)-(5)⇒ y ′ (xi ) = i+1 + R1 , kde |R1 | ≤ c1 h2 . 2h y − 2yi + yi−1 (4)+(5)⇒ y ′′ (xi ) = i+1 + R2 , kde |R2 | ≤ c2 h2 . h2 Konstanty c1 , c2 nezávisejí na h a i.
Máme tedy yi+1 − yi−1 + O(h2 ), 2h yi+1 − 2yi + yi−1 + O(h2 ). y ′′ (xi ) = 2 h y ′ (xi ) =
Dosadíme do a(x)y ′′ (x) + b(x)y ′ (x) − q(x)y(x) = f (x) v x = xi ; zavedeme ai = a(xi ), bi = b(xi ), atd. Dostaneme rovnost ai
yi+1 − 2yi + yi−1 yi+1 − yi−1 + bi − qi yi = fi + O(h2 ). 2h h2
Upravíme (stále jde o pˇresnou rovnost) −yi−1
ai − hbi /2 2ai + h2 qi ai + hbi /2 + y − yi+1 = −fi + O(h2 ). i 2 2 h h h2
Nyní zanedbáme cˇ len O(h2 ) a dostaneme sít’ovou rovnici pro pˇribližné uzlové hodnoty Yi−1 , Yi , Yi+1 : −Yi−1
ai − hbi /2 2ai + h2 qi a + hbi /2 + Y − Yi+1 i = −fi . i 2 2 h h h2
Rovnici sestavíme v každém vnitˇrním uzlu xi ∈ (a, b). Chyba diskretizace (rozdíl mezi −fi a levou stranou, v níž místo Yi−1 , Yi , Yi+1 užijeme hodnoty yi−1 , yi , yi+1 ) je ˇrádu O(h2 ). Jsou-li OP dirichletovské, pak hodnoty Y0 a YN známe a dostáváme N − 1 rovnic pro neznámé Y1 , . . . , YN−1 . Jsou-li v OP derivace, mužeme ˚ je aproximovat diferenˇcními podíly (Y1 − Y0 )/h a (YN − YN−1 )/h, dopustíme se však chyby ˇ diskretizace ˇrádu O(h), tedy vetší.
Chyba metody ηh = yh − Yh , kde yh = (y1 , . . . , yN )T a Yh = (Y1 , . . . , YN )T . Vektor Yh získáme vyˇrešením soustavy Ah Yh = Fh , kde Fh = (f1 , . . . , fN )T . Za pˇredpokladu hladkosti koeficientu˚ rovnice platí: ◮
matice Ah má "dobré" vlastnosti, existuje práveˇ jedno ˇrešení Yh ;
◮
dirichletovské OP, pak kηh kmax ≤ C0 h2
◮
◮
∀h ∈ (0, h0 ];
OP s derivací aproximované s chybou O(h), pak kηh kmax ≤ C1 h ∀h ∈ (0, h0 ];
ˇ OP s derivací aproximované s chybou O(h2 ) (tj. pˇresneji; 2 pˇredpis zde neuveden), pak kηh kmax ≤ C1 h ∀h ∈ (0, h0 ].
Poznámky: Chyba diskretizace OP je velmi významná. Snížení hladkosti ˇrešení a koeficientu˚ vede ke zhoršení odhadu˚ rychlosti konvergence.
ˇ Pˇríklad: Rešme y ′′ − (1 + sin2 x)y = −4, y(0) = 0 = y(π)
pˇribližneˇ metodou sítí s h = π/3. ˇ x0 = 0, x1 = π/3, x2 = 2π/3, √ Uzly síte: x3 = π. 2 2 q1 = 1 + sin x1 = 1 + sin (π/3) = 1 + ( √ 3/2)2 = 1 + 3/4, 2 2 q2 = 1 + sin x2 = 1 + sin (2π/3) = 1 + ( 3/2)2 = 1 + 3/4, f1 = −4, f2 = −4. Diferenˇcní rovnice Y0 − 2Y1 + Y2 π2 9
Y1 − 2Y2 + Y3 π2 9
− (1 +
3 )Y1 = −4, 4
(6)
− (1 +
3 )Y2 = −4. 4
(7)
Protože z okrajových podmínek plyne Y0 = 0 = Y3 , jsou (6)-(7) dveˇ rovnice pro dveˇ neznámé, tj. Y1 a Y2 .
Po odeˇctení −3(Y2 − Y1 )
7 (Y2 − Y1 ) = 0, 4 27 7 + (Y2 − Y1 ) = 0, 4 π2
π2 9
−
tedy Y1 = Y2 . Po dosazení do (6)-(7) Y1 = 1, 50
a
Y2 = 1, 50.
1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0
0
0.5
1
1.5
2
2.5
3
3.5
ˇ pˇribližné modˇre, 4 uzly. Pˇresné ˇrešení cˇ ervene,
1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0
0
0.5
1
1.5
2
2.5
3
3.5
ˇ pˇribližné modˇre, 5 uzlu. Pˇresné ˇrešení cˇ ervene, ˚
1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0
0
0.5
1
1.5
2
2.5
3
3.5
ˇ pˇribližné modˇre, 10 uzlu. Pˇresné ˇrešení cˇ ervene, ˚
1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0
0
0.5
1
1.5
2
2.5
3
3.5
ˇ pˇribližné modˇre, 20 uzlu. Pˇresné ˇrešení cˇ ervene, ˚
Pˇrednosti metody sítí: ◮
dosti obecná rovnice;
◮
jednoduchost; soustava s rˇídkou maticí.
◮
Zápory metody sítí: ◮
silné pˇredpoklady na hladkost;
◮
pˇri oslabení pˇredpokladu˚ na hladkost lze oˇcekávat pomalou konvergenci;
◮
problémy s okrajovými podmínkami (2D úlohy).
ˇ Problém vlastních císel a metoda sítí Víme, že vlastní cˇ ísla okrajové úlohy u ′′ + λu = 0, u(0) = 0,
u(π) = 0,
(8) (9)
ˇ jsou λ = 1, 4, 9, 16, 25, . . . Zkusme je spoˇcítat pˇribližne. ˇ x0 = 0, x1 = π/3, x2 = 2π/3, x3 = π, krok delení ˇ Uzly síte: h = π/3. Oznaˇcme U1 , U2 hodnoty sít’ového ˇrešení v bodech x1 , x2 a µ aproximaci vlastního cˇ ísla λ. Z (8) dostáváme diferenˇcní rovnice U0 − 2U1 + U2 + µU1 = 0, h2 U1 − 2U2 + U3 + µU2 = 0, h2 kde U0 = 0, U3 = 0, viz (9).
Po dosazení U0 = 0 a U3 = 0 soustavu −2U1 + U2 + µU1 = 0, h2 U1 − 2U2 + µU2 = 0, h2 upravíme 1 2 − 2 + µ U1 + 2 U2 = 0, h h 1 2 U1 + − 2 + µ U2 = 0. h2 h
2 − 2 +µ h
2
1 Soustava má nenulové ˇrešení, pokud − 4 = 0. h 2 1 1 3 Tedy µ − 2 = 2 , tj. µ1 = 2 ≈ 0, 9119 a µ2 = 2 ≈ 2, 7357. h h h h
ˇ Soustavu Jiný pˇrístup k soustave: −2U1 + U2 + µU1 = 0, h2 U1 − 2U2 + µU2 = 0, h2 mužeme ˚ upravit i takto 2 U1 − h2 1 − 2 U1 + h
1 U2 = µU1 , h2 2 U2 = µU2 ; h2
dostáváme (a pak ˇrešíme) standardní problém vlastních cˇ ísel 2 −1 “AU = µU” s maticí A = h12 . −1 2
ˇ delení: ˇ Jemnejší napˇr. h = π/5, vektor hodnot ve vnitˇrních uzlech síteˇ U = (U1 , . . . , U4 )T . Pak AU = µU, kde 2 −1 0 0 1 −1 2 −1 0 . A= 2 h 0 −1 2 −1 0 0 −1 2 Vlastní cˇ ísla (a vlastní vektory) tˇrídiagonální matice. ˇ ˇ delení, ˇ Cím jemnejší tím lepší aproximace (malých) vlastních cˇ ísel.
49
n=2 n=3 n=5 n = 10 n = 20 n = 50 n = 100
36 25 16 9 4 1 1
2 3 4 5 6 Poradove cislo vlastniho cisla
7
ˇ Je zobrazeno Parametr n udává poˇcet vnitˇrních uzlu˚ síte. ˇ nekolik malých vlastních cˇ ísel odpovídající matice typu n × n. S rostoucím n pozorujeme konvergenci k pˇresným vlastním cˇ íslum ˚ 1, 4, 9, 16, 25, 36, 49.
ˇ pˇrípad. Problém vlastních cˇ ísel a metoda sítí. Obecnejší ˇ Najdeme pˇribližneˇ první dveˇ vlastní cˇ ísla OÚ u ′′ + λp(x)u = 0, u(a) = 0,
u(b) = 0,
kde p je funkce kladná na intervalu [a, b]. ˇ x0 = a, x1 = (b − a)/3, x2 = 2(b − a)/3, x3 = b, krok Uzly síte: ˇ delení h = (b − a)/3. Zaved’me pi = p(xi ), i = 1, 2. Oznaˇcme U1 , U2 hodnoty sít’ového ˇrešení v bodech x1 , x2 a µ aproximaci vl. cˇ . λ. Diferenˇcní rovnice −2U1 + U2 + µp1 U1 = 0, h2 U1 − 2U2 + µp2 U2 = 0. h2
upravme
−2 1 + µp1 U1 + 2 U2 = 0, h2 h −2 1 U1 + + µp2 U2 = 0. h2 h2
Dveˇ hodnoty µ tedy najdeme jako ˇrešení kvadratické rovnice −2 −2 1 + µp1 + µp2 − 4 = 0. 2 2 h h h (Hodnoty p1, p2 a h jsou známé!)
Výchozí rovnice −2U1 + U2 + µp1 U1 = 0, h2 U1 − 2U2 + µp2 U2 = 0 h2 mužeme ˚ upravit i do této podoby 1 1 (2U1 − U2 ) = µU1 , p1 h 2 1 1 (−U1 + 2U2 ) = µU2 , p2 h 2 v níž již vidíme maticový vlastních cˇ ísel, tj. AU = µU, ! problém 1 0 2 −1 kde A = h12 p1 1 . −1 2 0 p 2
ˇ delení; ˇ ˇ ˇ Jemnejší maticový zápis. Stejnomerné delení s krokem h, body xi , v nich známé hodnoty pi a neznámé hodnoty Ui , i = 1, . . . , n. Pak AU = µU, kde U = (U1 , . . . , Un )T ,
a
2 −1 0 1 −1 A = 2D 0 h 0 0 0
p1 0 D= 0 0
−1 2 −1 ... ... ... 0
0 0 ... 0 −1 0 . . . 0 2 −1 0 . . . ... ... ... 0 −1 2 −1 0 0 −1 2 −1 . . . 0 −1 2
0 ... p2 . . . .. .. . . 0 0
0 0 . 0 pn
