Pokročilé numerické metody II
12. přednáška
Parabolické rovnice
Jiří Zelinka
Jiří Zelinka Pokročilé numerické metody II, 12. přednáška 1 / 14
Parabolické DR
Rovnice vedení tepla
u = u(x, t), f = f (x, t), x ∈ [0, 1], t ≥ 0
∂
∂t
u = α
∂2
∂x2
u + f (1)
Interpretace:
u(x, t) je teplota tenké tyče v daném bodě a daném čase,
funkce f zahrnuje vnější zdroje tepla,
α ≥ 0 je zpravidla konstanta – materiálové vlastnosti.
Počáteční podmínka: u(x, 0) = u0(x), x ∈ [0, 1]
Okrajové podmínky (Dirichletovy): u(i, t) = gi (t), t ≥ 0,
i ∈ {0, 1}
Homogenní rovnice: f ≡ 0.
Jiří Zelinka Pokročilé numerické metody II, 12. přednáška 2 / 14
Obecná rovnice
c(x, t)
∂u
∂t
+ q(x, t)u =
∂
∂x
(p(x, t)
∂u
∂x
) + f (x, t)
Newtonovy okrajové podmínky:
p(0, t)
∂u(0, t)
∂x
= α0u(0, t) − β0(t)
−p(1, t)
∂u(1, t)
∂x
= α1u(1, t) − β1(t)
Interpretace okrajových podmínek
u(0, t), u(1, t): teplota na krajích tyče, u Dirichletových
podmínek udržujeme stanovenou teloptu v daném čase.
∂u(0,t)
∂x
, ∂u(1,t)
∂x
: udává unikání tepla z konců tyče,
volný únik tepla: ∂u(0,t)
∂x
= 0.
Jiří Zelinka Pokročilé numerické metody II, 12. přednáška 3 / 14
Vlastnosti rovnice
Princip maxima
Každá funkce u(x, t), která splňuje homogenní rovnici
nabývá svého maxima i minima pro t = 0 nebo x ∈ {0, 1}.
Podmínky kompatibility (pro Dirichletovy podmínky)
u0(0) = g0(0), u0(1) = g1(0)
Rychlost šíření tepla
Pro x ∈ (−∞, ∞) (bez okrajových podmínek) se dá dokázat
u(x, t) =
1
2
√
απt
∞
−∞
u0(ξ)e−
(x−ξ)2
4αt dξ
Pro u0 ≥ 0 je funkce u je všude kladná, přestože u0 mohla být
kladná jen na omezeném intervalu ⇒ rychlost šíření tepla by
musela být nekonečná.
Jiří Zelinka Pokročilé numerické metody II, 12. přednáška 4 / 14
Zhlazovací jev
Pro spojitou a omezenou funkci u0 je řešení u dané
předchozím vztahem nekonečně hladké.
Jiří Zelinka Pokročilé numerické metody II, 12. přednáška 5 / 14
Přibližné řešení hledáme na [0, 1] × [0, T], T > 0.
Definice uzlů:
h = 1/N, x0 = 0, xN = 1, xi = i h,
τ = T/M, t0 = 0, tM = T, tk = k τ.
Označení: uk
i ≈ u(xi , tk).
Přibližné řešení počítáme postupně po časových vrstvách.
Diskretizaci rovnice lze provést různými způsoby.
Jiří Zelinka Pokročilé numerické metody II, 12. přednáška 6 / 14
Vyšetřování stability
Fourierova metoda: předpokládáme, že chybu v počáteční
podmínce můžeme vyjádřit pomocí komplexní Fourierovy
řady, šíření poruchy sledujeme pro jeden člen řady eiλx
.
Maticová analýza stability: zkoumáme vlastní čísla
transformačních matic. Pro α konstantní dostáváme tzv.
Toeplitzovy matice (matice s konstantními diagonálami),
které jsou navíc třídiagonální.
T =







d u 0 · · · 0
l d u · · · 0
... ... ...
0 · · · l d u
0 · · · 0 l d







typu n × n,
λk = d + 2
√
ul cos kπ
n+1
, k = 1, . . . , n.
Jiří Zelinka Pokročilé numerické metody II, 12. přednáška 7 / 14
A. Explicitní schéma
Rovnici aproximujeme v bodě (xi , tk):
1
τ
[uk+1
i − uk
i ] =
α
h2
[uk
i−1 − 2uk
i + uk
i+1] + f k
i
uk+1
i =
ατ
h2
[uk
i−1 − 2uk
i + uk
i+1] + uk
i + τf k
i
r =
ατ
h2
uk+1
i = r uk
i−1 + (1 − 2r)uk
i + r uk
i+1 + τf k
i ,
pro i = 1, . . . , N − 1, chyba aproximace je O(τ + h2
).
Pro r = 1
6
je chyba aproximace O(τ2
+ h4
).
Jiří Zelinka Pokročilé numerické metody II, 12. přednáška 8 / 14
Uk
=







uk
1
uk
2
...
uk
N−2
uk
N−1







, Fk
=







τf k
1 + ruk
0
τf k
2
...
f k
N−2
f k
N−1 + ruk
N







Uk+1
= Ar Uk
+ Fk
Ar =









1 − 2r r 0 0 · · · 0
r 1 − 2r r 0 · · · 0
0 r 1 − 2r r · · · 0
...
...
...
r 1 − 2r r
0 0 · · · 0 r 1 − 2r









Jiří Zelinka Pokročilé numerické metody II, 12. přednáška 9 / 14
Stabilita explicitního schématu
Schema je podmíněně stabilní. spektrální poloměr Ar musí být
menší nebo roven jedné, odtud
r ≤
1
2
t.j.
ατ
h2
≤
1
2
Důsledek
Časový krok závisí na druhé mocnině h.
Jiří Zelinka Pokročilé numerické metody II, 12. přednáška 10 / 14
B. Implicitní schéma
Rovnici aproximujeme v bodě (xi , tk+1):
1
τ
[uk+1
i − uk
i ] =
α
h2
[uk+1
i−1 − 2uk+1
i + uk+1
i+1 ] + f k+1
i
−r uk+1
i−1 + (1 + 2r)uk+1
i − r uk+1
i+1 = uk
i + τf k+1
i
Br Uk+1
= Uk
+ Fk+1
Br =









1 + 2r −r 0 0 · · · 0
−r 1 + 2r −r 0 · · · 0
0 −r 1 + 2r −r · · · 0
...
...
...
−r 1 + 2r −r
0 0 · · · 0 −r 1 + 2r









Schema je nepodmíněně stabilní, chyba je O(τ + h2
).
Jiří Zelinka Pokročilé numerické metody II, 12. přednáška 11 / 14
C. Crankovo–Nicholsonové schéma
Rovnici aproximujeme v bodě (xi , tk + τ/2):
1
τ
[uk+1
i − uk
i ] =
1
2
α
h2
[uk
i−1 − 2uk
i + uk
i+1]+
+
α
h2
[uk+1
i−1 − 2uk+1
i + uk+1
i+1 ] + f
k+1/2
i
−
r
2
uk+1
i−1 +(1+r)uk+1
i −
r
2
uk+1
i+1 =
r
2
uk
i−1+(1−r)uk
i +
r
2
uk
i+1+τf
k+1/2
i
Br/2Uk+1
= Ar/2 Uk
+ Fk+1/2
Schema je nepodmíněně stabilní, chyba je O(τ2
+ h2
).
Jiří Zelinka Pokročilé numerické metody II, 12. přednáška 12 / 14
Věta o konvergenci
Nechť řešení rovnice (1) s Dirichletovými okrajovými
podmínkami má řešení, které je čtyřikrát spojitě
diferencovatelné podle x a dvakrát podle t, 0 < r ≤ 1
2
. Pak
max |u(xi , tk) − uk
i | = O(τ + h2
)
pro explicitní a implicitní schéma a
max |u(xi , tk) − uk
i | = O(τ2
+ h2
)
pro schéma Crankovo–Nicholsonové.
Jiří Zelinka Pokročilé numerické metody II, 12. přednáška 13 / 14
Doplňky
Derivace v okrajových podmínkách:
Použijeme metodu fiktivních bodů podobně jako pro ODR.
Rovnice ve více dimenzích
∂
∂t
u = α∆u + f
Rovnici řešíme na Ω × R+, Ω ⊆ Rn
.
Např. pro n = 2 a explicitní schéma vyjde (viz eliptické
rovnice):
uk+1
ij = r uk
i,j−1 + r uk
i−1,j + (1 − 4r)uk
ij + r uk
i+1,j + r uk
i,j+1 + τf k
ij .
Podmínka stability: r ≤ 1
4
.
Jiří Zelinka Pokročilé numerické metody II, 12. přednáška 14 / 14