Pokročilé numerické metody I
Matice při numerickém řešení DR
Jiří Zelinka
Jiří Zelinka Matice při numerickém řešení DR 1 / 15
Obyčejné DR
Přibližné řešení hledáme v izolovaných bodech (uzlech)
x0, . . . , xN, zpravidla ekvidistantních, tj. xi+1 = xi + h.
Náhrada derivací (a chyby formulí)
u′
(xi ) =
1
h
[u(xi+1) − u(xi )] + O(h)
u′
(xi ) =
1
h
[u(xi ) − u(xi−1)] + O(h)
u′
(xi ) =
1
2h
[u(xi+1) − u(xi−1)] + O(h2
)
u′′
(xi ) =
1
h2
[u(xi−1) − 2u(xi ) + u(xi+1)] + O(h2
)
Odvození: z Taylorova rozvoje.
Jiří Zelinka Matice při numerickém řešení DR 2 / 15
Okrajová úloha: lineární rovnice 2. řádu
p(x)u′′
(x) + q(x)u′
(x) + r(x)u(x) = f (x), x ∈ [a, b]
Okrajové podmínky: u(a) = u0, u(b) = uN.
Definice uzlů: h = (b − a)/N, x0 = a, xN = b, xi = x0 + i h.
Označení: ui ≈ u(xi ), podobně pro p, q, r, f (přesné hodnoty).
Náhrada DR:
pi
1
h2
[ui−1 − 2ui + ui+1] + qi
1
2h
[ui+1 − ui−1] + ri ui = fi , tj.
(
pi
h2
−
qi
2h
)ui−1 + (ri −
2pi
h2
)ui + (
pi
h2
+
qi
2h
)ui+1 = fi
pro i = 1, . . . , N − 1,
u0 a uN jsu dány okrajovými podmínkami.
Výsledkem je systém lineárních rovnic s třídiagonální maticí.
Jiří Zelinka Matice při numerickém řešení DR 3 / 15
Příklad:
−u′′
(x) + u(x) = f (x), x ∈ [a, b], , u(a) = u0, u(b) = uN
Systém rovnic:
−
1
h2
[ui−1 − 2ui + ui+1] + ui = fi
t.j.
−ui−1 + (2 + h2
)ui − ui+1 = h2
fi , i = 1, . . . , N − 1
Dostali jsme systém lineárních rovnic s třídiagonální
diagonálně dominantní maticí.
Řešení: LU rozklad, Croutova metoda, pomocí inverze NE!
Jiří Zelinka Matice při numerickém řešení DR 4 / 15
2. Parciální DR: Rovnice vedení tepla
u = u(x, t), f = f (x, t), x ∈ [0, 1], t ≥ 0
∂
∂t
u = α
∂2
∂x2
u + f
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: u(i, t) = gi (t), t ≥ 0, i ∈ {0, 1}
Jiří Zelinka Matice při numerickém řešení DR 5 / 15
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 Matice při numerickém řešení DR 6 / 15
A. Explicitní schema
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
).
Jiří Zelinka Matice při numerickém řešení DR 7 / 15
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 Matice při numerickém řešení DR 8 / 15
Třídiagonální Toeplizovy matice
Tn =









α β 0
γ α β 0
0
...
...
...
...
...
...
...
... 0
0 γ α β
0 γ α









Vlastní čísla Tn: λj = α − 2
√
βγ cos jπ
n+1
, j = 1, . . . , n
Pro matici Ar : λj = 1 − 2r − 2r cos jπ
n+1
= 1 − 2r(1 + cos jπ
n+1
)
Aby |λj | ≤ 1 musí být r ≤ 1
2
.
Jiří Zelinka Matice při numerickém řešení DR 9 / 15
Stabilita explicitního schematu
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 Matice při numerickém řešení DR 10 / 15
B. Implicitní schema
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 Matice při numerickém řešení DR 11 / 15
C. Crankovo–Nicolsonové schema
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 Matice při numerickém řešení DR 12 / 15
3. Parciální DR: Eliptická rovnice
Poissonova rovnice na oblasti G v rovině:
∂2
u(x, y)
∂x2
+
∂2
u(x, y)
∂y2
= f (x, y)
spolu s okrajovými podmínkami na hranici G.
Oblast rozdělíme pravidelnou sítí bodů (xi , yj ) s krokem h v
obou směrech, ui,j ≈ u(xi , yj ).
Aproximace derivace:
∂2
u(x, y)
∂x2
+
∂2
u(x, y)
∂y2
=
1
h2
[u(x + h, y) + u(x − h, y)+
+u(x, y + h) + u(x, y − h) − 4u(x, y)] + O(h2
)
Jiří Zelinka Matice při numerickém řešení DR 13 / 15
Rovnice pro obdélníkovou oblast:
G = [a, b] × [c, d], b − a = N · h, d − c = M · h, a = x0,
c = y0
ui+1,j + ui−1,j + ui,j+1 + ui,j−1 − 4ui,j = h2
fi,j
i = 1, . . . , N − 1, j = 1, . . . , M − 1
Maticový zápis:
A · U = F
U =
(u1,1, . . . , uN−1,1, u1,2, . . . , uN−1,2, . . . , u1,M−1, . . . , uN−1,M−1),
podobně F s okrajovými podmínkami.
Jiří Zelinka Matice při numerickém řešení DR 14 / 15
A =









M I 0 0 · · · 0
I M I 0 · · · 0
0 I M I · · · 0
...
...
...
I M I
0 0 · · · 0 I M









,
I je jednotková matice řádu N − 1,
M =









−4 1 0 0 · · · 0
1 −4 1 0 · · · 0
0 1 −4 1 · · · 0
... ... ...
1 −4 1
0 0 · · · 0 1 −4









,
řád matice A je (N − 1) × (M − 1).
Jiří Zelinka Matice při numerickém řešení DR 15 / 15