ランダムウォーク 1
拡散現象
不可逆な過程
１ステップの移動
δ
θ
δ : 定数
θ: ランダム
拡散現象
ランダムウォーク 2
1000 粒子反射壁
１次元ランダムウォーク
−3δ −2δ −δ
0
(1) 各粒子は τ 秒ごとに距離
δ
2δ
3δ
δ だけ右か左に移動する.
(2) 各ステップで左右に行く確率はそれぞれ 1/2 であり,
前のステップでどちらに動いたかは記憶していない.
分布は空間的に一様化
注：個々の粒子は有意に空きスペースに行こうとする訳ではない.
(3) 各粒子は他の粒子と独立に動き, 相互作用することは
ない.
粒子の位置の平均
粒子の位置の分散
初期時刻に I 個の全粒子が原点にいる場合を考える.
第 i 番目の粒子の第 nステップにおける位置 :
xi (n)
−δ +δ
確率 1/2
xi (n) = xi (n − 1) ± δ
1
第 n ステップにおける粒子の平均位置 : !x(n)" = I
I
!
I
(xi (n
i=1
1
1) ± ) =
I
I
xi (n
i=1
1) = ⇥x(n
xi (n)
i=1
1)⇤
粒子の広がり具合
時間変数
t を導入する.
!
"
x(n)2 = nδ 2
δ2
D=
2τ
とおくと
標準偏差 σ(t)
σ(t) =
√
2Dt
x(n)
2
"
I
1#
=
xi (n)2
I i=1
"
xi (n)2 = xi (n
!x(n − 1)#
粒子数が十分大きいとき !x(n)" との関係は
1
⇤x(n)⌅ =
I
!
2
第 nステップにおける粒子の位置の分散 : x(n)
!
xi (n − 1)
!x(n)" = 0
!x(0)" = 0
1)2 ± 2 xi (n 1) + 2
!
"
!
" !
"
x(n)2 = x(n − 1)2 + δ 2 = x(n − 2)2 + 2δ 2
!
"
= · · · = x(0)2 + nδ 2
!
"
x(n)2 = nδ 2
粒子の広がる速さ
t = τn
!
" δ2
2
x(t) = t
τ
!
"
2
x(t) = 2Dt
粒子の広がり具合の指標
Ex. 4-1
空気中の微粒子の拡散係数
! 10−1 cm2 /sec
サイズ 10m 程度の部屋の
端から端まで香水の分子が
拡散によって移動するのに
どれぐらいかかるか？
D : 拡散係数
[D] = L2 T −1
長い距離を拡散で移動するには時間がかかる！
粒子数とゆらぎ
ランダムウォーク再び
1.0
Box 1
Box 2
粒子数：大
0.5
0.0
10粒子
t
密度ゆらぎ：小
1.0
0.5
0.0
100粒子
ランダムさを含まない
決定論的方程式で記述
可能になる．
t
1.0
0.5
比率
0.0
1000粒子
t
1.0
拡散のマクロな描像
0.5
0.0
拡散方程式（2 Box System）
Box 1
密度 :
u1 (t)
Box 2
密度 :
u1 (t)
u2 (t)
面積 : A
Ex.4-2
u2 (t)
面積 : A
Box 1
Box 2
単位時間あたり u1 (t) 移動
Box 1
Box 2
単位時間あたり u2 (t) 移動
u1 , u2 の時間発展を記述する方程式を求めよ．
du1
A
= (u2
dt
du2
u1 ), A
= (u1
dt
du1
= k(u2 − u1 )
dt
du2
= k(u1 − u2 )
dt
u2 )
k = /A とおいて
2 Box System
最も簡単な拡散方程式
10000粒子
t
2 Box System の方程式を解く
Ex.4-3
du1
= k(u2 − u1 )
dt
du2
= k(u1 − u2 )
dt
0
0
初期値を (u1 (0), u2 (0)) = (u1 , u2 ) とするとき
この方程式を解け.
Hint: ξ = u1 + u2 , η = u1 − u2
u2
とおいて ξ と η の方程式にする
d
d⇥
= 0,
=
dt
dt
(t)
(0) =
(t) =
0e
u1 (t) =
2k
u01
2kt
+ u02
= (u01
u02 )e
u01 + u02
u0 u02
+ 1
e
2
2
2kt
u1
2kt
, u2 (t) =
u01 + u02
2
u01
u1 + u2 が保存され, 最終的には一様になる
2
u02
e
2kt
,
拡散方程式（離散１次元版）
u1
ui−1 ui ui+1
du1
= k(u2 − u1 )
dt
dui
= k(ui−1 − ui ) + k(ui+1 − ui )
dt
duI
= k(uI−1 − uI )
dt
uI
離散ラプラシアン
dui
= k(ui−1 − ui ) + k(ui+1 − ui )
dt
= k(ui−1 − 2ui + ui+1 )
= 2k
i = 2, 3, · · · , I − 1
Ex.4-4
u1 + u2 + · · · + uI が保存されることを示せ.
!
ui−1 + ui+1
− ui
2
"
ui
離散ラプラシアン
隣接セルの平均値と自分自身の差を測っている.
隣接セルの平均値の方に行こうとするダイナミクス
空間的な凹凸をなくそうとするはず
両辺をすべて足しあわせれば明らか．
シミュレーション
離散から連続へ
拡散とは空間均一化のダイナミクスである.
I
最終的には
1!
lim ui (t) = u
¯0 =
ui (0)
t→∞
I i=1
ui+1
ui−1
拡散方程式（１次元連続版）
拡散方程式（１次元連続版）
Ex. 4-5
u
密度分布
J(a)
0
密度分布の勾配に比例
してものが流れる.
J(b)
a
!
b
∂u
流束 J = −D
∂x
x
Fick の法則
−J(b) + J(a) = D
Ex. 4-6
!
b
a
!
b
a
"
∂u
∂2u
−D 2
∂t
∂x
∂2u
∂u
=D 2
∂t
∂x
for 0 < x < !, t > 0
u(0, t) = u(!, t) = 0
for t > 0
u(x, 0) = f (x)
for 0 < x < !
境界条件
初期条件
B.C.：境界条件
u(x, t) ?
I.C.
B.C.
x
#
dx = 0
∂2u
∂u
=D 2
∂t
∂x
１次元拡散方程式
D ：拡散係数（必ず正）
初期値境界値問題 (ディリクレ条件)
B.C.
a
∂2u
dx であることを示せ.
∂x2
(a, b) は任意であったので
区間 udx = −J(b) + J(a)
t
b
u の満たすべき方程式を求めよ.
適当な区間 (a, b) でのものの出入りを考えると
d
dt
!
boundary condition
I.C.：初期条件
initial condition
初期値境界値問題 (ノイマン条件)
∂2u
∂u
=D 2
∂t
∂x
for 0 < x < !, t > 0
∂u
∂u
(0, t) =
(", t) = 0
∂x
∂x
for t > 0
境界条件
u(x, 0) = f (x)
for 0 < x < !
初期条件
ノイマン境界条件は, 境界において物の出入り
がないことを意味している.
Ex.4-7
u
d
dt
! !
udx
の総量は保存されることを示せ.
Z
`
udx =
0
Z
0
`
0
@u
dx =
@t
Z
`
0

@2u
@u
D 2 dx = D
@x
@x
`
=0
x=0
初期値が正弦波の解
境界条件について
Ex.4-8
! ディリクレ (Dirichlet) 条件
u(0, t) = u(!, t) = 0 のように境界で
u の値が与えられる.
! ノイマン (Neumann) 条件
∂u
∂u
(0, t) =
(", t) = 0
∂x
∂x
のように境界で
初期値が正弦波 sin
(ディリクレ条件）
um (x, t) = am (t) sin
∂u
の値が与えられる.
∂x
! 周期境界条件これは実は境界がないという条件である.
∂u
∂u
! ! u(0,! t) = u(!,
! t), ∂x
(0, t) = ∂x (!, t)
円周上の関数
初期値が正弦波の解
Ex.4-9
πmx
初期値が正弦波 cos
であるような解を構成せよ.
"
(ノイマン条件）
πmx
(m = 0, 1, 2, · · ·)
um (x, t) = am (t) cos
"
! πm "2
dam
am
= −D
dt
l
πm
"
2
)
am (t) = e−D(
πm
"
2
πm
l
um (x, t) = e−D(
(m = 1, 2, 3, · · ·)
2
)
πm
"
t
) t sin πmx
"
2
モード数の２乗に比例して, 素早く減衰する.
シミュレーション
ディリクレ条件
4モード解
t
) t cos πmx
"
モード数の２乗に比例して, 素早く減衰する.
ただし...
um (x, t) = e−D(
⇡mx
`
! πm "2
dam
am
= −D
dt
l
am (0) = 1
am (0) = 1
am (t) = e−D(
πmx
であるような解を構成せよ.
"
8モード解
線形重ね合わせ
線形重ね合わせ
区間 [0, !] 上で拡散方程式の初期値境界値問題を考える.
u1 (x, t), u2 (x, t) がそれぞれ初期条件
u1 (x, 0) = f1 (x), u2 (x, 0) = f2 (x) を満たす解である時
それらの１次結合
u = a 1 u 1 + a2 u 2
u(x, 0) = f (x)
ただし f (x) = a1 f1 (x) + a2 f2 (x)
は初期条件
を満たす初期値境界値問題の解である．
Ex.4-11
∂2u
∂u
= D 2 の初期値
区間 (0, ) 上で, 拡散方程式 ∂t
∂x
境界値問題を考える. ただし, 境界条件はディリクレ
x
2 x
u(x, 0) = 2 sin
sin のとき
条件で初期値が ⇥
⇥
u(x, t) を求めよ.
u1 (x, t) = exp
Ex.4-10
上のことを確かめよ.
D 2t
⇥2
sin
x
, u2 (x, t) = exp
⇥
と置くと、解は u1 と u2 の一次結合
2
2
u = 2u1
sin
2 x
⇥
u2 となる.
2
u2
u
u1
u2
u1
u
= a1
+ a2
= a1 D
+ a2 D
=D 2
t
t
t
x2
x2
x
u(0, t) = a1 u1 (0, t) + a2 u2 (0, t) = 0
4D 2 t
⇥2
u( , t) = a1 u1 ( , t) + a2 u2 ( , t) = 0
u(x, t) = 2 exp
D 2t
⇥2
sin
x
⇥
exp
4D 2 t
⇥2
sin
2 x
⇥
u(x, 0) = a1 u1 (x, 0) + a2 u2 (x, 0) = a1 f1 (x) + a2 f2 (x) = f (x)
線形重ね合わせ
Ex.4-12
∂2u
∂u
=
D
区間 (0, ) 上で, 拡散方程式 ∂t
∂x2 の初期値
境界値問題を考える. ただし, 境界条件はノイマン条件
1
2 x
3 x
u(x, 0) =
cos
+ 2 cos
で初期値がのとき
3
⇥
⇥
u(x, t) を求めよ.
u0 (x, t)
1
u2 (x, t) = exp
4D 2 t
⇥2
cos
2 x
, u3 (x, t) = exp
⇥
と置くと、解は定数と u2 とu3 の一次結合 u =
1
u(x, t) =
3
exp
4D 2 t
⇥2
2 x
cos
+ 2 exp
⇥
シミュレーション
1
u0
3
9D 2 t
⇥2
cos
3 x
⇥
u2 + 2u3 となる.
9D 2 t
⇥2
3 x
cos
⇥
4モード解
＋
20モード解
まず、高いモードの解（短波長の解）が先に
減衰して、低いモードの解はゆっくりと .....
一般の初期値について (ディリクレ条件)
f (x)
初期値をフーリエサイン展開して構成する.
f (x) =
∞
!
bm sin
m=1
Ex.4-13
初期値境界値問題の解
πmx
"
u(x, t) を級数の形で求め，
lim u(x, t) = 0
初期値が何であれであることを示せ.
t→∞
u(x, t) =
D
bm exp
2
m2 t
sin
⇥2
m=1
mx
⇥
フーリエ級数を用いて解く
Ex.4-14
∂2u
∂u
=D 2 の
区間 (0,1) 上で, 拡散方程式 ∂t
∂x
初期値境界値問題を考える. ただし, 境界条件は
ディリクレ条件で初期値は次のように与えられる.
u(x, 0) = x(1
フーリエ級数を用いて
x(1
x) =
8
3
n:odd
(0
x)
u(x, t) を求めよ.
1
sin nx
n3
8 X 1
u(x, t) = 3
e
⇡
n3
lim u(x, t) = 0 は明らか．
t→∞
1)
x
D⇡ 2 n2 t
(0
x
1)
sin ⇡nx
n:odd
一般の初期値について (ノイマン条件)
f (x)
初期値をフーリエコサイン展開して構成する.
f (x) = a0 +
∞
!
am cos
m=1
Ex.4-15
初期値境界値問題の解
πmx
"
u(x, t) を級数の形で求め，
lim u(x, t) を求めよ.
t→∞
u(x, t) = a0 +
am e
m=1
lim u(x, t) = a0 =
t→∞
1
!
!
0
D(
m
⇥
!
f (x)dx
2
) t cos mx
⇥
f (x) の平均値
フーリエ級数を用いて解く
Ex.4-16
∂2u
∂u
=D 2 の
区間 (0,1) 上で, 拡散方程式 ∂t
∂x
初期値境界値問題を考える. ただし, 境界条件は
ノイマン条件で初期値は次のように与えられる.
!
1 for 0 < x < 1/2
u(x, 0) =
0 for 1/2 < x < 1
フーリエ級数を用いて u(x, t) を求めよ.
また lim u(x, t) は何になるか.
t→∞
u(x, 0) =
1
2
+
2
1
2
u(x, t) = +
2
n:odd
n:odd
sin ( n/2)
cos nx
n
sin ( n/2)
e
n
D
2
n2 t
cos nx
lim u(x, t) =
t
⇥
1
2
数値解法：差分法
偏微分方程式の数値解法の中で, 元の
微分方程式の微分の部分を差分で置き
換えて解く方法を差分法と言う.
拡散方程式の
シミュレーション
他にも，有限要素法，有限体積法，境界要素法，
スペクトル法など，いろいろある．
U " (x) の差分近似
Ex.4-17
U(x) を x のまわりで３次の項までTaylor 展開せよ.
関数 (0 < δx ! 1)
U !! (x) 2 U !!! (x) 3
U ! (x)
δx +
δx +
δx + O(δx4 )
1!
2!
3!
€
U !! (x) 2 U !!! (x) 3
U ! (x)
U (x − δx) = U (x) −
δx +
δx −
δx + O(δx4 )
1!
2!
3!
U (x + δx) = U (x) +
€
Ex.4-18
(1) U ! (x) を U (x + δx) と U (x) を用いて近似せよ.
U " (x) =
U(x + δx) − U(x)
+ Ο(δx)
δx
U ! (x) !
前進差分, 右側差分
€
y
U ! (x) !
U " (x) =
U(x) − U(x − δx)
+ Ο(δx)
δx
(3) U ! (x) を U (x + δx)と U (x − δx) を用いて近似せよ.
€
U(x + δx) − U(x − δx)
U " (x) =
+ Ο(δx 2 )
2δx
y = U (x)
U (x) − U (x − δx)
δx
後退差分, 左側差分
(2) U ! (x) を U (x − δx) と U (x) を用いて近似せよ.
€
U (x + δx) − U (x)
δx
U ! (x) !
U (x + δx) − U (x − δx)
2δx
中心差分
x − δx
x
x + δx
x
U !! (x) の差分近似
区間の離散化のしかた
Ex.4-19
U !! (x) を U (x), U (x + δx), U (x − δx) を用いて近似せよ.
U !! (x) =
U (x − δx) − 2U (x) + U (x + δx)
+ O(δx2 )
δx2
U (x − δx) − 2U (x) + U (x + δx)
U (x) !
δx2
京都方式
xi = iδx (i = 0, 1, · · · , I)
δt
tn = nδt
(n = 0, 1, 2, · · ·)
xi = (i − 1/2)δx
xi
(i = 1, 2, · · · , I)
x
空間：札幌, 時間：京都
xi = (i − 1/2)δx (i = 1, 2, · · · , I)
時間微分は前進差分で，
空間微分は中心２階差分で
2
u
u
=D 2
t
x
δx
tn
札幌方式
拡散方程式の差分化
空間と時間の離散化
t
!
δx = "/I
!!
中心２階差分, 離散ラプラシアン
I 等分
0
uni = u(xi , tn )
u(x, t + t)
t
u(x, t)
D
u(x
x, t)
uni−1 − 2uni + uni+1
un+1
− uni
i
=D
δt
δx2
2u(x, t) + u(x + x, t)
x2
t
(i = 1, 2, · · · , I; n = 0, 1, · · ·)
x
境界条件の処理
境界条件の処理
n
n
差分のcellを両側に1つずつ広げて，そこでの値を u0 , uI+1 とおく.
unI+1
un0
! ノイマン条件
un1
unI
un
u
(0, tn ) = 0, I+1
x
x
un0
x
u
( , tn ) = 0
x
un0 = un1 , unI+1 = unI
unI
un1
! 周期境界条件
! ディリクレ条件
un0 + un1
2
u(0, tn ) = 0,
unI + unI+1
2
0 1 2
u( , tn ) = 0
I 1 I I+1
n
n
n
un
0 = −u1 , uI+1 = −uI
n
n
n
un
0 = uI , uI+1 = u1
初期値境界値問題の計算スキーム
拡散方程式
un − 2uni + uni+1
un+1
− uni
i
= D i−1
δt
δx2
(i = 1, 2, · · · , I; n = 0, 1, · · ·)
境界条件
ディリクレ条件
ノイマン条件
周期境界条件
un
0
un0
un
0
=
=
=
n
n
−un
1 , uI+1 = −uI
un1 , unI+1 = unI
n
n
un
I , uI+1 = u1
初期条件
u0i = f (xi )
(i = 1, 2, · · · , I)
新しいステップの計算
uni−1 − 2uni + uni+1
un+1
− uni
i
=D
δt
δx2
λ=
Dδt
とおくと
δx2
un+1
= uni + λ(uni−1 − 2uni + uni+1 )
i
Explicit scheme
(n = 0, 1, · · ·)
n+1
n
i−1
i
i+1
計算手順
von Neumann の安定性解析
計算の初め
I.C.
u0i = f (xi ) (i = 1, 2, · · · , I)
B.C.
u00 , u0I+1
D.Eq.
数値シミュレーションでは, 各計算ステップで微小な
誤差が混入する. これが計算ステップとともに指数的
に成長するようならば, 計算は破綻する！
波数
u1i (i = 1, 2, · · · , I)
√
−1kx の成長を見る.
k の周期関数 e
境界条件は無視.
√
−1kxi
uni = µn e
成長係数
µ
とおいて差分方程式に代入し,
を調べる.
すべての波数
k に対し |µ| ≤ 1
ある波数 k に対し |µ| > 1
Explixcit Scheme の安定条件
Ex.4-20
拡散方程式の Explicit scheme の安定性を調べよ.
√
µn+1 e
−1kxi
√
− µn e
δt
−1kxi
=D
Dδt
(2 cos kδx − 2)
µ−1=
δx2
k は任意の数
安定条件
δt ≤
δx2
2D
1−
√
µn e
−1kxi−1
√
√
− 2µn e −1kxi + µn e
δx2
4Dδt
≤µ≤1
δx2
−1 ≤ 1 −
4Dδt
δx2
空間の解像度を上げるために δx を
小さくとると, δt を非常に小さく
取らなければならない.
不安定
安定条件が破れると
un+1
i
=
uni
−1kxi+1
kδx
4Dδt
sin2
µ=1−
δx2
2
安定
2Dδt
+
δx2
2Dδt
>1
δx2
un+1
i
!
uni−1 + uni+1
− uni
2
"
un+1
i
uni−1
uni−1 + uni+1
は
2
を行き過ぎる！
uni
un+1
i−1
uni+1
un+1
i+1
メッシュごとのジグザグ波形で指数発散
これが出たら, 安定性の破れをまず疑うべし！
Fully Implicit Scheme
レポート課題
un − 2uni + uni+1
uni − un−1
i
= D i−1
δt
δx2
n
後退差分
Dδt
λ=
とおくと
δx2
Ex.
uni
n
1 + (1 + 2 )ui
上式を行列で Aun = un
Dirichlet
1+2
A=
..
..
..
..
i−1
uni+1 = uni
..
..
..
A=
..
1+
..
..
..
1+3
Fully Implicit scheme
un − 2uni + uni+1
uni − un−1
i
= D i−1
δt
δx2
Crank-Nicolson scheme
uni − un−1
D
i
=
δt
2
..
..
..
..
..
..
1+2
1+2
1
と書くとき，を書き下せ。
A
1+2
..
i+1
i
1
Neumann
1+3
次の２つの拡散方程式の計算スキームの安定性を調べよ.
1
n
1+
!
n−1
un−1
+ un−1
uni−1 − 2uni + uni+1
i−1 − 2ui
i+1
+
δx2
δx2
"