1 次元拡散方程式

1 次元拡散方程式

11 次元拡散方程式

拡散方程式の数値解法
一次元拡散方程式∂U∂ 2U=κ 2∂t∂x(1)について考える.

(1) 式で境界条件として U の値を与える場合をディリクレ問題,場合をノイマン問題という.

∂Uの値を与える∂x空間方向に離散化して考える. xn = n∆x とおいて 0 ≤ x ≤ L の領域を N 個に離散化する (図 1 参照).

L
∆x0x1 x2 x3xnxN

図 1: 0 ≤ x ≤ L の領域を N 個に離散化する. xn = n∆x とする.

Lである. (1) 式の U を上記の離散表現と三角関数を用いて波数ごNとに離散化すると,∑U (xn , t) =Ũl (t)e−ikl (xn )ここで, ∆x =l=∑Ũl (t)e−ikn (n∆x)(2)lとなる. (2) 式を (1) 式に代入すると∑ d Ũll

2012 0202-ogihara.tex

dte−ikl (n∆x) = −κ∑kl2 Ũl (t)e−ikl (n∆x) .

l

2013/06/05(荻原 弘尭)

1 次元拡散方程式

1 次元拡散方程式2

両辺に e−ikn (n∆x) を掛けて 0 から 2π まで積分を取ると,∫ 2π ∑0ld Ũl −ikl (n∆x) −ikn (n∆x)
eedx = −κdt

∫ 2π ∑0kl2 Ũl (t)e−ikl (n∆x) e−ikn (n∆x) dx.

l

三角関数の直行性から n 成分以外が 0 になるので,

d Ũn (t)= −κkn2 Ũn (t)dt(3)となる. ここで周期境界条件を仮定すると kn はkn =2πn2πn=N ∆xL(4)となる. 波を表現できる最大の波数 kn は n = N のときなので,kmax =2π∆x(5)となる. (3) 式は κkn2 = K と置くと, 摩擦係数が K の摩擦方程式と同じ形になる.
ただし, 安定条件を考える場合には波数によって安定条件が変わってくるので最も波数の大きいモードについて考えなければならない.

各スキームの安定性条件
反復しない 1 段階スキームを用いた場合と, 2 段階スキームの代表であるリープフロッグスキームの場合を考える.

(3) 式で kn = kmax としたとする. その式を時間方向に 1 段階スキームを用いて離散化すると,(Ũn+1= Ũ − κn2π∆x)2∆t(αŨ n+1 + β Ũ n ),(6)α+β =1となる. ここで時間差分スキーム (2) の (1) 式を用いた. 摩擦方程式で考察したとき( 2π )2と同様に K = κ ∆x∆t とおいて安定性を議論する.

1. α = 1, β = 0(オイラースキーム) の場合

2012 0202-ogihara.tex

2013/06/05(荻原 弘尭)

1 次元拡散方程式

1 次元拡散方程式

3

摩擦方程式からオイラースキームの安定性は0<K≤1となる. よって,

(0<κ2π∆x)2∆t ≤ 1(7)2. α = 0, β = 1(後退スキーム) の場合
摩擦方程式より,

1≤11+K(8)となり, 全ての K で安定である. よって, ∆t, ∆x によらず安定である. しかしながら, 実際にこのスキームで解く場合 K が 1 程度の値じゃないと良い精度の解は出ない.

3. α = β = 12 (台形スキーム) の場合
摩擦方程式より,0<K≤2となる. よって,

(0<κ2π∆x)2∆t ≤ 2(9)4. リープフロッグスキームの場合(3) 式で kn = kmax としたとする. その式をリープフロッグスキームを用いて離散化すると,()22πn+1n−1Ũ
= Ũ
− 2κ∆tŨ n(10)∆x( 2π )2となる. 1 段階スキームと同様に K = κ ∆x∆t とおいた場合, 摩擦方程式と同様の形になるので摩擦方程式と同じ安定性条件になる.
リープフロッグスキームの摩擦方程式における安定性を考える. 摩擦方程式
からリープフロッグスキームの増幅率は√λ1 = −K + K 2 + 1√
λ2 = −K − K 2 + 1
2012 0202-ogihara.tex

2013/06/05(荻原 弘尭)

1 次元拡散方程式1 次元拡散方程式4である. ここでそれぞれ物理モードと計算モードになる. 物理モードの場合は常に安定性条件を満たす. しかし, 計算モードは常に不安定でありステップ毎に符号が変わる.
解の中から計算モードを完全に取り除くことができない以上リープフロッグスキームは摩擦方程式に向かない.
安定性条件は摩擦方程式と同じだが ∆x の二乗にもよることに注意が必要である.

問題例
実際に初期条件と境界条件を与えて一次元拡散方程式を考えていく. 考える方程式系は以下である.

∂U∂ 2U= κ 2,∂t∂x0 ≤ x ≤ 1,(11)初期条件 : U (x, 0) = sin (πx),
境界条件 : U (0, t) = U (1, t) = 0.

解析解(11) 式の解析解を変数分離法によって求める. U (x, t) をU (x, t) = X(x)T (t)(12)と置く. (12) 式を (11) 式に代入すると,X(x)dT
d2X= κT (t) 2 .
dtdx

両辺を X(x)T (t) で割ると,

1 d2X
1 dT
=.
κT (t) dt
X(x) dx2

2012 0202-ogihara.tex

2013/06/05(荻原 弘尭)

1 次元拡散方程式

1 次元拡散方程式5

左辺は t だけの関数, 右辺は x だけの関数となっている. x と t は互いに独立である
ので左辺と右辺が等号で結ばれるにはこの値が定数にならなければならない. 分離定数を a とすると,

1 dT= κaT (t) dt
1 d2X=aX(x) dx2(13)(14)となる. (14) 式から

d2X
= aX
dx2

(15)

X(x) = An cos(kn x) + Bn sin(kn x)(16)となる. 今とすると, 境界条件よりX(0) = An = 0,X(1) = Bn sin(kn ) = 0.
よって,kn = nπ(n = 0, 1, 2, · · · ).

よって, (16) 式はX(x) = Bn sin(nπ)(17)となる. (15) 式と同じになるには

a = −(nπ)2(18)となる. (18) 式を (13) 式に代入すると,

dT= −κ(nπ)2 T.
dt
これを解くと,T (t) = T (0)e−κ(nπ) t .
2(19)(19) 式と (17) 式より一般解はU (x, t) =∞∑Cn sin(nπx)e−κ(nπ) t2(20)n=02012 0202-ogihara.tex

2013/06/05(荻原 弘尭)

1 次元拡散方程式

1 次元拡散方程式6となる. ここで T (0) と Bn はともに定数なので一つにまとめて Cn とした. 初期条件 U (x, 0) = sin(πx) より,U (x, 0) =∞∑Cn sin(nπx)

n=0= sin πx.
この式が成り立つにはC1 = 1,C0 = 0,
Cn>1 = 0.
よって, 解析解は

U (x, t) = sin(πx)e−κ(π) t2(21)となる.

数値解法(11) 式を数値解法で解く. 時間方向にはオイラー法, 空間方向には中心差分を用いて離散化する.
解U を

U (i∆x, n∆t) = Uinと置く. このとき (11) 式を時間方向にはオイラー法, 空間方向には中心差分を用いて離散化する. 2 階微分が出ているので半格子を図 2 のように入れて考える.

L∆x0x1 x2 x3∆xxi−1 xi

xi+1xNxi− 1 xi+ 122

図 2: 0 ≤ x ≤ L の領域を N 個に離散化する. i を挟むように格子点の半分の値を半格子として定義する.

2012 0202-ogihara.tex

2013/06/05(荻原 弘尭)

1 次元拡散方程式

1 次元拡散方程式

7

このとき,

Uin+1 − Uin1=κ∆t∆x(∂U n∂U n−∂x i+ 1∂x i− 1)22( n)
nUi+1 − Uin Uin − Ui−11=κ
−∆x
∆x
∆xnnn
U − 2Ui + Ui−1= κ i+1.
(∆x)2よって, Uin+1 について解くと,)∆t ( nnn
U−2U+Ui+1ii−1(∆x)2(22)U0n = UNn = 0.

(23)

Uin+1 = Uin + κとなる. ただし, 境界条件より実際の数値解の振る舞い(22) 式を実際に数値的に解いてみた. 図 3 に初期条件の状態を示す.

図 3: t = 0 における U (x, 0) の分布. 点線は U (x, 0)e−1 の分布を示す.

2012 0202-ogihara.tex

2013/06/05(荻原 弘尭)

1 次元拡散方程式

1 次元拡散方程式8

図 4: ∆x = 0.1, ∆t = 1, t = 10 での図. 実線は数値解, 破線は厳密解, 点線はU (x, 0)e−1 の分布を示す.

実線は数値解, 破線は厳密解, 点線は e-folding time での値を示している. 図 4 は∆x = 0.1, ∆t = 1 として数値的に t = 10 まで計算した図である.
減衰している様子がわかる. 今回の e-folding time は1κπ 2≈ 10τ=である. 図 4 の値は若干減衰が大きいがほぼ厳密解に等しい.
安定条件は境界条件によって異なる. 今回は固定端である. 固定端でのオイラー法の安定条件は( π )20 < κ∆t≤1∆xである. よって, 今回はκ∆t( π )2∆x= 0.01 × 1 × π 2 × 10−2
≈ 0.99となるので安定条件を満たしている.
2012 0202-ogihara.tex

2013/06/05(荻原 弘尭)

1 次元拡散方程式

1 次元拡散方程式

9K が (∆x)2 に反比例しているのでむやみやたらに ∆x を小さくとると ∆t をかなり小さくとらなければならなくなる. よって, 計算時間を減らすためには ∆x はあまり細かくとることはできない.

2012 0202-ogihara.tex

2013/06/05(荻原 弘尭)

1 次元拡散方程式

1 次元拡散方程式10

参考文献
石岡 圭一, 2004,「スペクトル法による数値計算入門」, 東京大学出版会, ISBN:4130613057
荻原 弘尭, 2010,「スペクトル法を用いた数値計算– 一次元線形移流方程式の場合–」
URL:http://www.ep.sci.hokudai.ac.jp/∼psg/doc2011/ogihara B/ogihara B.pdf
荻原 弘尭, 2011, 「時間差分スキーム (2)」
URL:http://www.ep.sci.hokudai.ac.jp/g̃fdlab/comptech/resume/0804/2011 0803-ogihara.pdf
荻原 弘尭, 2012, 「摩擦方程式」
URL:http://www.ep.sci.hokudai.ac.jp/g̃fdlab/comptech/resume/0112/2012 0112-ogihara.pdf

2012 0202-ogihara.tex

2013/06/05(荻原 弘尭)