lZ
0
[pu′
v′
+ quv]dx +
X
i
αiu(ci)v(ci) =
lZ
0
fvdx +
X
i
βiv(ci)
a(u, v) = L(v) ∀v
−(pu′
)′
+ qu = f
x0, . . . xN , wi(xj) = δijv
U(x) u(x), U(x) =
NX
i=0
Uiwi(x)
±p(ci)u′
(ci) = αiu(ci) − βi, ci ∈ {0, l}u(ci) = gi
aproximace řešení -> Metoda konečných prvků
Rovnice:
+ okr. podmínky nebo
Slabá formulace:
Volba funkcí : báze prostoru lin. splajnů na uzlech
Ik
(·)
[xk−1, xk]
lZ
0
pu′
v′
dx =
NX
k=1
xkZ
xk−1
pu′
v′
dx ≈
NX
k=1
Ik
(pu′
v′
),
lZ
0
quvdx ≈
NX
k=1
Ik
(quv)
Ik
(pu′
w′
k) = pk−1/2
Uk − Uk−1
hk
1
hk
hk = pk−1/2
Uk − Uk−1
hk
wk
[xk−1, xk] wk−1Obdélníkové pravidlo na pro funkci :
Ik
(pu′
w′
k−1) = pk−1/2
Uk − Uk−1
hk
(−
1
hk
)hk = pk−1/2
Uk−1 − Uk
hk
jsou aproximace integrálů přes jednotlivé subintervaly obdélníkovým nebo
lichoběžníkovým pravidlem.
Obdélníkové pravidlo na pro funkci :
[xk−1, xk] wkLichoběžníkové pravidlo na pro funkci :
[xk−1, xk] wk−1Lichoběžníkové pravidlo na pro funkci :
Ik
(pu′
w′
k) =
pk−1 + pk
2
Uk − Uk−1
hk
1
hk
hk =
pk−1 + pk
2
Uk − Uk−1
hk
Ik
(pu′
w′
k−1) =
pk−1 + pk
2
Uk−1 − Uk
hk
Podobně pro další integrály
Obdélníkové pravidlo:
Ik
(quwk) = qk−1/2Uk−1/2
1
2
hk = qk−1/2
Uk − Uk−1
4
hk
Ik
(quwk−1) = qk−1/2Uk−1/2
1
2
hk = qk−1/2
Uk − Uk−1
4
hk
Ik
(fwk) = fk−1/2
1
2
hk = Ik
(fwk−1)
Lichoběžníkové pravidlo:
Ik
(quwk−1) =
1
2
[qk−1Uk−11 + qkUk0]hk =
1
2
qk−1Uk−1hk
Ik
(quwk) =
1
2
[qk−1Uk−10 + qkUk1]hk =
1
2
qkUkhk
Ik
(fwk) =
1
2
[fk−10 + fk1]hk =
1
2
fkhk
Ik
(fwk−1) =
1
2
[fk−11 + fk0]hk =
1
2
fk−1hk
Ik
(pu′
w′
j)Doporučení - pro použít obdélníkové pravidlo, jinak lichoběžníkové.
Pro každou funkci získáme jednu lineární rovnici.wk
Ik
(pu′
w′
k) + Ik+1
(pu′
w′
k) + Ik
(quwk) + Ik+1
(quwk) + αkUk =
= Ik
(fwk) + Ik+1
(fwk) + βk
pro k = 1, . . . , N − 1
k = 0 p1/2
U0 − U1
h1
+
1
2
h1q0U0 + α0U0 =
1
2
f0h1 + β0
k = N pN−1/2
UN − UN−1
hN
+
1
2
hN qN UN + αN UN =
1
2
fN hN + βN
pro :
pro :
pk−1/2
Uk − Uk−1
hk
+ pk+1/2
Uk − Uk+1
hk+1
+
1
2
Uk[hkqk + hk+1qk] + αkUk =
=
1
2
fk[hk + hk+1] + βk
Matice soustavy pro neznámé U0, . . . , UN :
\frac{p_{1/2}}{h_1}+\frac{1}{2}h_1q_0+\alpha_0
p1/2
h1
+
1
2
h1q0 + α0 −
p1/2
h1
−
p1/2
h1
−
p3/2
h2
−
p3/2
h2
0
0 0
0
−
p5/2
h3
0
0
0
. . .
. . .
. . .
...
pk−1/2
hk
+
pk+1/2
hk+1
+
1
2
[hkqk + hk+1qk+1] + αkPrvky na hlavní diagonále:
Vedlejší prvky: −
pk−1/2
hk
, −
pk+1/2
hk+1
p1/2
h1
+
p3/2
h2
+
1
2
[h1q1 + h2q1] + α1
Diferenční metoda:
−(pu′
)′
+ qu = f
+ Dirichletovy okrajové podmínky, hledáme přibližné řešení na uzlech
x0, . . . xN
u′
(xk) ≈
1
2h
[u(xk+1) − u(xk−1)]
Náhrada derivací (pro ekvidistantní uzly):
u′′
(xk) ≈
1
h2
[u(xk+1) − 2u(xk) + u(xk−1)
−pk
1
h2
[uk+1 − 2uk + uk−1] − p′
k
1
2h
[uk+1 − uk−1] + qkuk = fk
Náhrada rovnice, pokud neznáme p':
zk+1/2 = pk+1/2u′
k+1/2 ≈ pk+1/2
uk+1 − uk
h
, zk−1/2 = pk−1/2u′
k−1/2 ≈ pk−1/2
uk − uk−1
h
Celkem
−
pk+1/2uk+1 − (pk+1/2 + pk−1/2)uk + pk−1/2uk−1
h2
+ qkuk = fk
uj xjNáhrada rovnice, pokud známe p' ( značí přibližnou hodnotu řešení v ):
k = 1, . . . , N − 1, u0 a uNpro jsou známé.
z = pu′
, (zk)′
≈
zk+1/2−zk−1/2
h
h2
Po vynásobení dostáváme rovnice
−pk+1/2uk+1 + (pk+1/2 + pk−1/2)uk − pk−1/2uk−1 + h2
qkuk = h2
fk
u0 a uN
p3/2 + p1/2 + h2
q1 −p3/2
−p3/2
p5/2 + p3/2 + h2
q2
−p5/2
−p5/2
−p7/2p7/2 + p5/2 + h2
q3
0 0
0
0
0
0
0
. . .
...
s maticí soustavy (členy s převedeme na pravou stranu)
a pravou stranou
[h2
f1 + p1/2u0, h2
f2, h2
f3, . . . , h2
fN−2, h2
fN−1 + pN−1/2uN ]T
hk = h ∀k
Porovnání s rovnicemi slabé formulace:
Pro Dirichletovy okrajové podmínky budou chybět diskrétní části
a dále předpokládáme ekvidistantní uzly, tj. . Rovnice
se tedy zjednoduší na
pk−1/2
Uk − Uk−1
hk
+ pk+1/2
Uk − Uk+1
hk+1
+
1
2
Uk[hkqk + hk+1qk] + αkUk =
=
1
2
fk[hk + hk+1] + βk
a po vynásobení dostanemeh
pk−1/2[Uk − Uk−1] + pk+1/2[Uk − Uk+1] + Ukh2
qk = fkh2
pk−1/2
Uk − Uk−1
h
+ pk+1/2
Uk − Uk+1
h
+
1
2
Uk[hqk + hqk] =
1
2
fk[h + h]
což po úpravě dá
k = 1, . . . , N − 1pro , tedy stejný systém rovnic jako pro diferenční metodu.
−pk−1/2Uk−1 + (pk−1/2 + pk+1/2)Uk − pk+1/2Uk+1 + h2
qkUk = h2
fk
Domácí úkol:
Program (v Matlabu) pro řešení rovnice -(pu')'+qu=f s Dirichletovými okrajovými
podmínkami.
Vstup: p, q, f, I=[a,b], N (pro ekvidist. dělení intervalu) nebo uzly
Výstup: U_k - přibližné hodnoty řešení v uzlech
Poslat do pondělí.