第11章時間発展問題の数値解法

はじめに
時間積分法
非定常熱伝導問題
θ法
まとめ
参考文献
.
.
第 11 章時間発展問題の数値解法
畔上秀幸
名古屋大学情報科学研究科複雑系科学専攻
July 17, 2014
.
.
.
.
.
.
1 / 19
はじめに
時間積分法
非定常熱伝導問題
θ法
まとめ
参考文献
§11.1 はじめに
(目標) 時間発展問題に対する数値解法について考えてみよう．ま
ず，常微分代数方程式の初期値問題に対する時間積分法を理解してか
ら，非定常の熱伝導問題に対する解法について考える．
.
.
.
.
.
.
2 / 19
はじめに
時間積分法
非定常熱伝導問題
θ法
まとめ
参考文献
§11.2 時間積分法
次のような常微分代数方程式の初期値問題を考える．
.
問題 11.2.1 (常微分代数方程式の初期値問題)
.
ψ : Rn → Rn が与えられたとき，
dy
= ψ (y) ,
dt
y (0) = y0
n
.を満たす y : (0, tT ) → R を求めよ．
ここでは， Runge-Kutta 法を用いて数値積分を行うことを考えよう．
.
.
.
.
.
.
3 / 19
はじめに
時間積分法
非定常熱伝導問題
θ法
まとめ
参考文献
§11.2 時間積分法 (cnt.)
.
定義 11.2.2 (2 次の Runge-Kutta 法)
.
問題 11.2.1 に対して，h > 0 を定数として， tk+1 = tk + h, t0 = 0,
k = 1, 2, · · · , とする．このとき，y (t) の近似解を
y (tk+1 ) = ψ (y (tk )) +
k1 + k2
h
2
で計算する方法を 2 次の Runge-Kutta 法という．ただし，
k1 = ψ (y (tk )) ,
k2 = ψ (y (tk ) + k1 h)
とする．
.
.
.
.
.
.
.
4 / 19
はじめに
時間積分法
非定常熱伝導問題
θ法
まとめ
参考文献
§11.2 時間積分法 (cnt.)
2 次の Runge-Kutta 法の誤差は h2 のオーダーである．なぜならば，
y (tk+1 ) = y (tk ) + y˙ (tk ) h +
( )
1
y¨ (tk ) h2 + O h3
2
に対して， a, b を未知数として
y (tk+1 ) = y (tk ) + aψ (y (tk )) h + bψ (y (tk ) + y˙ (tk ) h) h
= y (tk ) + aψ (y (tk )) h
)
(
( 2 ))
∂ψ (
h
ψ (y (tk )) y˙ (tk ) h + O h
+ b ψ (y (tk )) +
∂y T
= y (tk ) + aψ (y (tk )) h
(
(
( )))
h
+ b ψ (y (tk )) + ψ (y (tk )) + y¨ (tk ) h + O h2
( )
= y (tk ) + (a + b) ψ (y (tk )) h + bψ (y (tk )) + O h3
が成り立つことから，a = b = 1/2 となるためである．
.
.
.
.
.
.
5 / 19
はじめに
時間積分法
非定常熱伝導問題
θ法
まとめ
参考文献
§11.2 時間積分法 (cnt.)
.
定義 11.2.3 (4 次の Runge-Kutta 法)
.
問題 11.2.1 に対して，h > 0 を定数として， tk+1 = tk + h, t0 = 0,
k = 1, 2, · · · , とする．このとき，y (t) の近似解を
y (tk+1 ) = ψ (y (tk )) +
k1 + 2k2 + 2k3 + k4
h
6
で計算する方法を 4 次の Runge-Kutta 法という．ただし，
)
(
h
,
k1 = ψ (y (tk )) , k2 = ψ y (tk ) + k1
2
(
)
h
k3 = ψ y (tk ) + k2
, k4 = ψ (y (tk ) + k3 h)
2
とする．
.
.
.
.
.
.
.
6 / 19
はじめに
時間積分法
非定常熱伝導問題
θ法
まとめ
参考文献
§11.3 非定常熱伝導問題
問題 1.1 を再記しよう．ただし，cv = 1 および λ = I (d 次の単位行
列) とおく．
.
問題 11.3.1 (熱伝導問題)
.
b0 : (0, t¯) × Ω → R, pN : (0, t¯) × (∂Ω \ ΓD ) → R, uD : (0, t¯) × ΓD → R,
u0 : Ω → R が与えられたとき，
∂u
− ∆u = b0 in (0, t¯) × Ω,
∂t
∂ν u = pN on (0, t¯) × ΓN ,
u = uD on (0, t¯) × ΓD ,
u = u0 in Ω at t = 0
を満たす
u : (0, t¯) × Ω → R を求めよ．
.
.
.
.
.
.
.
7 / 19
はじめに
時間積分法
非定常熱伝導問題
θ法
まとめ
参考文献
§11.4 θ 法
τ を時間刻みを表す正定数とする．時間領域 (0, t¯) を n = t¯/τ に分割
し，k ∈ {0, 1, · · · , n} に対して tk = kτ とかく．uh0 ∈ Uh を u0 の近似
関数，k ∈ {0, 1, · · · , n} に対して bhk と phk を b0 (tk ) と pN (tk ) の近
似値とする．ある θ ∈ [0, 1] を選び，
( · )h k+θ = θ( · )h k+1 + (1 − θ) ( · )hk
のようにかく．このとき，
1
(uh k+1 − uhk , vh ) + a (uh k+θ , vh ) = lh k+θ (vh )
τ
(11.4.1)
を満たすように uh1 , · · · , uhn を求める方法を θ 法という．ただし，
∫
(uh , vh ) =
uh vh dx
Ωh
∫
∫
lh k+θ (vh ) =
bh k+θ vh dx +
ph k+θ vh dγ
Ωh
とする．
ΓNh
.
.
.
.
.
.
8 / 19
はじめに
時間積分法
非定常熱伝導問題
θ法
まとめ
参考文献
§11.4 θ 法 (cnt.)
θ = 0 のとき，前進 Euler 法，θ = 1 のとき，後退 Euler 法，θ = 1/2
のとき，Crank-Nicolson 法という．
空間領域 Ω に対する有限要素分割と近似関数の作り方については第
5 章に示したとおりとする．すなわち，k ∈ {0, 1, · · · , n} に対して，
(
) (
)
∑
∑
¯ kD
u
ϕkD
¯ k · ϕ,
¯ =
ukj ϕj =
·
=u
ukj ϕj +
uhk (u)
¯ kN
u
ϕkN
vhk (¯
v) =
j∈ND
j∈NN
∑
∑
j∈ND
vkj ϕj +
j∈NN
(
vkj ϕj =
(11.4.2)
) ( )
v¯kD
ϕD
·
= v¯k · ϕ
v¯kN
ϕN
(11.4.3)
とおく．
.
.
.
.
.
.
9 / 19
はじめに
時間積分法
非定常熱伝導問題
θ法
まとめ
参考文献
§11.4 θ 法 (cnt.)
(11.4.2) と (11.4.3) を (11.4.1) に代入すれば，
{
}
)
(
)
(
¯ − (1 − θ) τ A
¯ u
¯ + θτ A
¯ u
¯ k + τ θl¯k+1 + (1 − θ) l¯k
¯ k+1 = M
M
(11.4.4)
¯ は第 5 章で定義した全体係数行列，l¯k は第 5 章
とかける．ただし，A
で定義した全体既知項ベクトルにおいて b0 と pN をそれぞれ bhk と
phk におきかえたときのベクトルとする．また，
∑
¯ i Zi ∈ R|N |×|N | ,
¯ =
M
ZiT M
i∈E
)
(∫
¯i =
M
∈ R|Ni |×|Ni |
φi(α) φi(β) dx
Ωi
αβ
である．
.
.
.
.
.
.
10 / 19
はじめに
時間積分法
非定常熱伝導問題
θ法
まとめ
参考文献
§11.4 θ 法 (cnt.)
θ = 0 (前進 Euler 法) のとき，(11.4.4) は，
¯u
¯ −1 A
¯ k + τ l¯k
¯ k+1 = u
¯k − τ M
u
(11.4.5)
¯ −1 を計算しておけば，t0 = 0 のときの u
¯ 0 と l¯0
となる．このとき，M
¯ 1 が計算される．同様に k を更
を与えれば，逆行列計算をしなくても u
¯2, · · · , u
¯ n が計算される．このように，tk のデータ
新していくことで u
から tk+1 のデータを求めていく計算方法を陽解法という．
.
.
.
.
.
.
11 / 19
はじめに
時間積分法
非定常熱伝導問題
θ法
まとめ
参考文献
§11.4 θ 法 (cnt.)
それに対して，θ = 1 (後退 Euler 法) のとき，(11.4.4) は，
}
) {
(
¯u
¯ + τA
¯ −1 M
¯ k + τ l¯k+1
¯ k+1 = M
u
(11.4.6)
となる．このように，tk と tk+1 の既知データから tk+1 の未知データ
を求めていく計算方法を陰解法という．
.
.
.
.
.
.
12 / 19
はじめに
時間積分法
非定常熱伝導問題
θ法
まとめ
参考文献
§11.4 θ 法 (cnt.)
さらに，θ = 1/2 (Crank-Nicolson 法) のとき，
(
)−1 {(
)
)}
τ (¯
¯ + τA
¯ − τA
¯
¯ u
¯ k+1 = M
¯k +
M
lk+1 + l¯k
u
2
2
2
(11.4.7)
となる．
.
.
.
.
.
.
13 / 19
はじめに
時間積分法
非定常熱伝導問題
θ法
まとめ
参考文献
§11.4 θ 法 (cnt.)
θ 法の安定性について，次のことが知られている [1]．
Vh を要素の最大長さが h となる有限要素分割における近似関数の集
合とする．ある正定数 c0 が存在して，任意の vh ∈ Vh に対して
c0
∥vh ∥L2 (Ωh ;R)
∥vh ∥H 1 (Ωh ;R) ≤
h
が成り立つとき，Vh は逆不等式を満たすという．ただし，u : Ω → R を
領域 Ω ⊂ Rd 上の関数とする．p ∈ [1, ∞] および s ∈ [0, ∞] に対して，
∑ ∫ ∇β u (x)p dx < ∞ (1 ≤ p < ∞) ,
|β|≤s
Ω
|β|≤s
x∈Ω
max ess sup ∇β u (x) < ∞
(p = ∞)
を満たすような u の全体集合を Sobolev 空間といい，W s,p (Ω; R) とか
T
d
く．β = {β1 , · · · , βd } ∈ {0, · · · , d} は
∇β u =
∂ β1 ∂ β2 · · · ∂ βd u
∂xβ1 1 ∂xβ2 2
· · · ∂xβd d
,
|β| =
∑
i∈{1,··· ,d}
.
βi ≤ s
.
.
.
.
.
14 / 19
はじめに
時間積分法
非定常熱伝導問題
θ法
まとめ
参考文献
§11.4 θ 法 (cnt.)
を意味する多重指数である．W s,p (Ω; R) は
|u|W i,p (Ω;R)

1/p

∑



∇β u (x)p p

(1 ≤ p < ∞)
L (Ω;R)
=
|β|=i



max ∇β u (x) ∞
(p = ∞)
L (Ω;R)
|β|=i
をセミノルムとして
∥u∥W s,p (Ω;R) =
(
s
∑


p

|u|
)1/p
W i,p (Ω;R)


 max |u(x)| i,∞
W
(Ω;R)
(1 ≤ p < ∞)
i=0
0≤i≤m
(p = ∞)
をノルムとして実 Banach 空間になる．W 0,2 (Ω; R) を L2 (Ω; R) とも
かく．W 1,2 (Ω; R) を H 1 (Ω; R) ともかく．
.
.
.
.
.
.
15 / 19
はじめに
時間積分法
非定常熱伝導問題
θ法
まとめ
参考文献
§11.4 θ 法 (cnt.)
.
定理 11.4.1 (θ 法の安定性)
.
Vh は逆不等式を満たすとする．θ ∈ [0, 1/2) のときは，ある δ ∈ (0, 1)
に対して
τ≤
2(1 − δ) 2
h
(1 − 2θ)c20
が成り立つように時間刻み τ をとるとする．このとき，θ 法の解 uh0 ,
· · · , uhn に対して




∑
1
2
2
¯lk+θ 2 1′
∥uhk ∥L2 (Ωh ;R) ≤
∥uh0 ∥L2 (Ωh ;R) + τ
H (Ωh ;R) 
δ
k∈{0,··· ,n−1}
が成り立つ．ただし，
H 1′ (Ωh ; R) は H 1 (Ωh ; R) の双対空間である．
.
.
.
.
.
.
.
16 / 19
はじめに
時間積分法
非定常熱伝導問題
θ法
まとめ
参考文献
§11.4 θ 法 (cnt.)
.
θ ∈ [1/2, 1] のときには τ を任意にとる．このとき，θ 法の解 uh0 , · · · ,
uhn に対して
∑
2
2
¯lk+θ 2 1′
|uhk | ≤ ∥uh0 ∥L2 (Ω;R) + τ
H (Ω;R)
k∈{0,··· ,n−1}
.が成り立つ．
.
.
.
.
.
.
17 / 19
はじめに
時間積分法
非定常熱伝導問題
θ法
まとめ
参考文献
§11.5 まとめ
非定常の熱伝導問題は θ 法により解けることをみた．
.
.
.
.
.
.
18 / 19
はじめに
時間積分法
非定常熱伝導問題
θ法
まとめ
参考文献
参考文献
[1] 田端正久.
偏微分方程式の数値解析.
岩波書店, 東京, 2010.
.
.
.
.
.
.
19 / 19

Download Report