はじめに 時間積分法 非定常熱伝導問題 θ法 まとめ 参考文献 . . 第 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
© Copyright 2025 ExpyDoc