u

講義予定(案)
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
( 9/ 2) 数値シミュレーションの手続き (テキスト第１章)
( 9/ 9) 偏微分方程式と解析解 (テキスト第２章)
( 9/16) 休講
( 9/30) 差分方程式とそのスキーム (テキスト第３章) +変換 (テキスト第４章)
(10/ 7) 計算 (テキスト第５章)+連立一次方程式の解法(テキスト第６章)
(10/21) 流れ関数‐ポテンシャルによる解法(テキスト第７章)
(10/28) 流速‐圧力を用いた解法 (テキスト第７章)
(11/ 4) 熱流体解析と多相流解析
(11/11) 乱流の数値解析 by 金子暁子先生
(11/18) 数値解析の実際 by 渡辺正先生(JAEA)
(11/25) (予備日)
変換ー変化する方程式の定性的性質ー
工学への応用
(開発、設計、性能評価)
変換
(差分近似)
(離散化)
仮説
(定式化)
物理現象
(熱流体挙動)
数学モデル
(偏微分方程式)
数値計算
計算モデル
(差分方程式)
適合性
適切性
安定性
Laxの同等定理
厳密解
収束性
差分解
代数方程式
数値解
有限差分方程式とそのスキーム
空間(x,y)と時・空間(x,t)の離散化
Δx
Δx
終端値
j+1
n+1
Δy
j
j-1
y
Δt
n
n-1
t
i-1
境界点
i
i+1
内点
x
0
i-1
i
i+1
初期値
x
差分法における空間の記述
Δx
連続的な空間(x,y)を等間隔に
x方向にΔx
y方向にΔy
j+1
で分割。
j
x方向にi番目、
y方向にj番目
おける変数を
u(iΔx,jΔy)
j-1
y
ui , j
と記述する。
Δy
i-1
境界点
i
i+1
内点
x
差分法における時空間の記述
連続的な時・空間(x,t)を等間隔に
Δx
空間をΔx、
終端値
時間をΔt
で分割。
x方向にi番目、
n番目の時間
における変数
u(iΔx,nΔt)を
n
i
u
と記述する。
n+1
Δt
n
n-1
t
0
i-1
i
i+1
初期値
x
空間の前進差分近似（1/2）
ui+1,jを点(i,j)においてxについてTaylor展開し
ui +1, j
3
1 ∂ 2u
1
u
∂u
∂
3
2
x
x
= ui , j +
Δx +
Δ
+
Δ
+L
3
2
2 ∂x i , j
6 ∂x i , j
∂x i , j
∂u / ∂x i , j について解く
ui + 1 , j − ui , j 1 ∂ 2 u
∂u
2
Δ
x
+
O
Δ
x
−
=
Δx
2 ∂x 2 i , j
∂x i , j
( )
=
ui + 1 , j − ui , j
Δx
+ O (Δx )
空間の前進差分近似（2/2）
これは、導関数 ∂u / ∂x を表現する差分商 Δu / Δx の1つが
∂u Δu ui +1, j − ui , j
≈
=
+ O (Δx )
∂x Δ x
Δx
であることを示している。これを空間の前進差分近似
（forward space difference）という。
Taylor展開の打ち切りの状況から、この式は空間について
「1次の正確度を持つ」という。
空間の後退差分近似（1/2）
ui-1,jを点(i,j)においてxについてTaylor展開
ui −1, j
3
∂
∂u
1 ∂ 2u
1
u
2
3
Δ
−
Δ
+L
= ui , j −
Δx +
x
x
2
3
∂x i , j
2 ∂x i , j
6 ∂x i , j
∂u / ∂x i , j について解く
∂u
∂x
=
i, j
=
ui , j − ui −1, j
Δx
ui , j − ui −1, j
Δx
1 ∂ 2u
+
2 ∂x 2
+ O (Δx )
( )
Δx + O Δx 2
i, j
空間の後退差分近似（2/2）
これは、導関数 ∂u / ∂x を表現する差分商 Δu / Δx の1つが
∂u Δu ui , j − ui − 1 , j
≈
=
+ O (Δx )
∂x Δ x
Δx
であることを示している。これを空間の後退差分近似
（backward space difference）という。
Taylor展開の打ち切りの状況から、この式は空間について
1次の正確度を持つ。
空間の中心差分近似（1/2）
ui+1,j 、ui-1,jを点(i,j)においてxについてTaylor展開
∂u
1 ∂ 2u
1 ∂ 3u
2
3
= ui , j +
Δx +
Δ
+
Δ
+L
x
x
2
3
∂x i , j
2 ∂x i , j
6 ∂x i , j
ui +1, j
ui −1, j
∂u
1 ∂ 2u
1 ∂ 3u
2
= ui , j −
Δx +
Δx −
Δx 3 + L
2
3
∂x i , j
2 ∂x i , j
6 ∂x i , j
上式から下式を差し引く
ui +1, j − ui −1, j
1 ∂ 3u
∂u
3
Δ
+L
x
Δx + 2 ⋅
= 2⋅
6 ∂x 3 i , j
∂x i , j
∂u / ∂x i , j について解く
∂u
∂x
=
i, j
ui + 1, j − ui −1, j
2Δx
( )
+ O Δx 2
空間の中心差分近似（2/2）
これは、導関数 ∂u / ∂x を表現する差分商 Δu / Δx の1つが
( )
∂u Δu ui + 1, j − ui −1, j
≈
=
+ O Δx 2
∂x Δ x
2Δx
であることを示している。これを空間の中心差分近似
（centered space difference）という。
Taylor展開の打ち切りの状況から、この式は空間について
2次の正確度を持つ。
２階の導関数の差分商（1/2）
ui+1,j 、ui-1,jを点(i,j)においてxについてTaylor展開
ui +1, j
∂u
1 ∂ 2u
1 ∂ 3u
2
= ui , j +
Δx +
Δx +
Δx 3 + L
2
3
∂x i , j
2 ∂x i , j
6 ∂x i , j
ui −1, j
∂u
1 ∂ 2u
1 ∂ 3u
2
3
= ui , j −
Δx +
Δ
−
Δ
+L
x
x
2
3
∂x i , j
2 ∂x i , j
6 ∂x i , j
上式と下式を足し合わせる
ui +1, j + ui −1, j
∂ 2 u / ∂x 2
∂ 2u
∂x
i, j
について解く
i,j
=
2
∂ 2u
= 2ui , j + 2 Δx 2 + L
∂x i , j
ui + 1, j − 2ui , j + ui −1, j
Δx
2
( )
+ O Δx 2
２階の導関数の差分商（2/2）
これは、導関数 ∂u / ∂x を表現する差分商 Δu / Δx の1つが
( )
Δu ui + 1, j − 2ui , j + ui −1, j
2
O
x
≈
=
+
Δ
∂x 2 Δ x
Δx 2
∂ 2u
であることを示している。これを空間の中心差分近似
（centered space difference）という。
Taylor展開の打ち切りの状況から、この式は空間について
2次の正確度を持つ。
空間の差分商（1/2)
ui , j の時空間の点(i,ｊ)における空間についての差分商
∂u uin+1 − uin
=
+ O (Δx )
Δx
∂x
（前進差分近似）
∂u uin − uin−1
+ O (Δx )
=
Δx
∂x
（後退差分近似）
( )
∂u uin+ 1 − uin−1
=
+ O Δx 2
2Δx
∂x
∂ 2u
∂x
2
=
uin+ 1 − 2uin + uin−1
Δx
2
（中心差分近似）
( )
+ O Δx 2
時空間の差分商（2/2)
u in の時空間の点(i,n)における時間についての差分商
∂u uin+1 − uin
=
+ O (Δt )
∂t
Δt
（前進差分近似）
∂u uin − uin −1
=
+ O (Δt )
∂t
Δt
（後退差分近似）
( )
∂u uin + 1 − uin −1
=
+ O Δt 2
2 Δt
∂t
∂ 2u
∂t
2
=
uin + 1 − 2uin + uin −1
Δt
2
（中心差分近似）
( )
+ O Δt 2
波動方程式(双曲型)の差分スキーム
波動方程式の有限差分スキーム
波動方程式
2
∂ 2u
2 ∂ u
−c
=0
2
2
∂t
∂x
差分方程式(中心差分)
n
n
n
uin +1 − 2uin + uin −1
2 ui + 1 − 2ui + ui −1
2
2
−
c
+
O
Δ
t
+
O
Δ
x
=0
2
2
Δt
Δx
( ) ( )
uin+1 について解く。
Δt
λ =c
u = − u + λ u + 2(1 − λ )u + λ u
Δx
→ CFL（Courant-Friedrichs-Lewy）スキーム
λ：CFL数（CFL number）
n+1
i
n −1
i
2
n
i +1
2
n
i
2
n
i −1
波動方程式の差分スキームのプログラム例
uin+1についての差分方程式
(
)
uin+1 = − uin −1 + λ2 uin+1 + 2 1 − λ2 uin + λ2 uin−1
Δt
λ =c
Δx
uin+1 を求めるためのプログラム (BASIC)
1490
1500
1510
1520
*CFL
FOR I=2 TO IM1
U(I)=-U2(I)+L2*U1(I+1)+2*(1-L2)*U1(I)+L2*U1(I-1)
NEXT I:RETURN
対流方程式の（陽）差分スキーム
対流方程式：
∂u
∂u
=0
+c
∂x
∂t
差分方程式：
uin+1 − uin
uin+1 − uin−1
+c
+ O (Δt ) + O Δx 2 = 0
Δt
2Δx
( )
uin+1 について解く。
uin + 1
=
uin
(
C n
−
ui + 1 − uin−1
2
)
Δt
C=c
Δx
→ FTCS（Forward-Time Centered Space）スキーム
C：Courant数（Courant number）
クーラン(Courant)条件

速度をcとすると、微小時間⊿t に進む距離は、c ・⊿tと
なる。物理的には、この距離が、差分した微小区間⊿xを
超えると計算が発散してしまう。そのためには、
cΔt < Δx
でなければならない。すなわち、以下の条件を、Courant
条件といい、計算結果の安定性に影響を与えものとであ
る。
cΔt
<1
cΔt
Δx
Δx
対流方程式の差分スキーム(その１)
（FTCSスキーム）
Δx
uin+1 についての差分方程式
uin + 1 = uin −
(
C n
ui + 1 − uin−1
2
)
Δt
C=c
Δx
n+1
Δt
n
t
n-1
FTCSスキームによるプログラム例
i-1
1620 *FTCS
1630 FOR I=2 TO IM1
1640 TP1=-.5*CDTDX*(UN(I+1)-UN(I-1))
1660 U(I)=UN(I)+TP1:NEXT I:RETURN
i
i+1
x
対流方程式の差分スキーム(その２)
蛙とび法（Leap frog scheme）
Δx
蛙とび法による差分方程式
uin + 1
=
uin −1
−C
(
uin+ 1
− uin−1
)
n+1
Δt
n
n-1
t 0
蛙とび法によるプログラム例
1670 *LPFRG
1680 FOR I=2 TO IM1
1690 TP1=-CDTDX*(UN(I+1)-UN(I-1))
1710 U(I)=U2(I)+TP1:NEXT I:RETURN
i-1
i i+1
x
対流方程式の差分スキーム(その3)
風上差分（Upwind scheme）
Δx
風上差分法による差分方程式
(
− C (u
uin + 1 = uin − C uin − uin− 1
n+ 1
i
u
=u
n
i
n
i+1
−u
n
i
)
)
C≥0
Δt
n
C≤0
風上差分法によるプログラム例
1570
1580
1590
1610
n+1
n-1
t
0
*UPWND
FOR I=2 TO IM1
TP1=-CDTDX*(UN(I)-UN(I-1))
U(I)=UN(I)+TP1:NEXT I:RETURN
i-1
i i+1
C≥0
x
対流方程式の陰解法(Implicit scheme) (1/3)

対流方程式：
∂u
∂u
=0
+c
∂t
∂x
対流方程式を時刻n+1において、時間：後退差分近似、
空間：中心差分近似すると、次のような差分方程式が得
られる。
uin+ 1 − uin
uin++11 − uin−+11
= −c
Δt
2 Δx
n+ 1
uin++11 、 uin + 1 、 ui − 1 を変数値として差分スキームを求
めると、
1 n+ 1
1 n+ 1
n+ 1
αui + 1 + ui − αui − 1 = uin
2
2
対流方程式の陰解法(Implicit scheme) (2/3)

0 ≤ x ≤ 1 の区間をi分割した場合を考えると： t
i=1のとき
−

i=2のとき
−

α
2
α
i=3のとき
−

i=iのとき
−

2
u0n + 1 + u1n + 1 +
u1n + 1 + u2n + 1 +
α
2
α
2
α
2
α
u2n + 1 + u3n + 1 +
uin−+11 + uin + 1 +
2
2
n+1
n
u3n + 1 = u2n
2
α
Δt
u2n + 1 = u1n
α
u4n + 1 = u3n
uin++11 = −
以上の連立方程式を行列で表すと、
Δx
α
0
1
…
α⎞
⎛
uin−+11 + ⎜ 1 + ⎟ uin + 1 = uin
2
2⎠
⎝
i-1
i
x
対流方程式の陰解法(Implicit scheme) (3/3)

差分方程式の行列式：
α
⎡
⎢ 1
⎢ α
⎢−
⎢ 2
⎢ 0
⎢
⎢
⎢ M
⎢
⎢ M
⎢
⎢
⎢ 0
⎣

2
1
−
α
2
O
0
α
L
2
O
1
O
O
O
O −
L
L
α
2
0
L
O
α
2
1
−
α
2
⎤
0⎥
⎡ n α n+ 1 ⎤
⎥ n+ 1
u + u
M ⎥ ⎡ u1 ⎤ ⎢ 1 2 0 ⎥
⎥
n
⎥ ⎢ u n+ 1 ⎥ ⎢
u
2
2
⎢
⎥ ⎢
⎥
M ⎥⎥ ⎢ M ⎥ ⎢
⎥
M
⎥
⎥⎢ M ⎥ = ⎢
M
⎥ ⎢
⎥
0 ⎥⎢
n
1
n
+
⎥
⎥ ⎢ ui − 2 ⎥ ⎢
ui − 1
α ⎥ ⎢ n+ 1 ⎥ ⎢
⎥
α
n
n+ 1
u
⎢
⎥
2 ⎥ ⎣ i − 1 ⎦ ⎢ ui − ui ⎥
2
⎣
⎦
⎥
1⎥
⎦
境界条件
∂uin
=0
∂x
を利用して決定する
以上の行列の式を解くことによって、各時刻における変位が求まる。
Fi+1/2
最新の空間差分スキーム
- QUICK -
i+1
セルiで積分
１次元対流方程式
∂φ
∂ (uφ )
=−
∂t
i
φ
∂x
n+ 1
i
= φ − Δt
n
i
Fi n+ 1 2 − Fi n− 1 2
Vi
Fi n+ 1 2 = uin+ 1 2φin+ 1 2 S i + 1 2
セル境界の物理量 uin−+11/ 2 を境界の前後を挟む２つのセルと
風上側に１つ離れたセルの値を用いて２次の多項式で内挿
⎧1
⎪ (− φ i − 1 + 6φ i + 3φ i + 1 ) u ≥ 0
⎪8
⎪ 1 (− φ + 6φ + 3φ ) u < 0
i+2
i+1
i
⎪⎩ 8
φi + 1 2 = ⎨
差分式：
φ
n+ 1
i
3φi + 1 + 3φi − 7 φi − 1 + φi − 2
⎧ n
tu
φ
Δ
−
⎪⎪ i
8 Δx
=⎨
⎪φ n + Δtu 3φi − 1 + 3φi − 7 φi + 1 + φi + 2
⎪⎩ i
8 Δx
u≥0
u<0
最新の空間差分スキーム
- QUICKEST -
Fi+1/2
i
i+1
uin−+11/ 2
セル境界の物理量
を境界の前後を挟む２つのセルと
風上側に２つ離れたセルの値を用いて２次の多項式で内挿
2
c
c
n
n
n
n
⎧ φ − (2φ i +1 + 3φ i − 6φ i −1 + φ i − 2 ) + (φ in+1 − 2φ in + φ in−1 )
⎪
6
2
⎪
c3 n
⎪
− (φ i +1 − 3φ in + 3φ in−1 − φ in− 2 ) c ≥ 0
6
⎪
=⎨
c
c2 n
n
n
n
n
n
⎪ φ i + (2φ i −1 + 3φ i − 6φ i +1 + φ i + 2 ) + (φ i −1 − 2φ in + φ in+1 )
6
2
⎪
3
⎪
c
+ (φ in−1 − 3φ in + 3φ in+1 − φ in+ 2 ) c < 0
⎪
6
⎩
n
i
φ in +1
差分式：
φi + 1 2
⎧1
⎪⎪ (φ i + φ i + 1 ) −
2
=⎨
⎪ 1 (φ + φ ) −
i+1
⎪⎩ 2 i
c
1 − c2
(φi + 1 − φi ) −
(φi + 1 − 2φi + φi − 1 )
2
6
c
1 − c2
(φi + 1 − φi ) −
(φi + 2 − 2φi + 1 + φi )
2
6
最新の空間差分スキーム
- TVD (Total Variation Diminishing) (1/2)セル境界で物理量を内挿する際、数値振動が発生する場所には１次精度の
風上差分を、それ以外では高次精度の差分スキームを用いる。
φin+ 1 2
(
)
⎧ n 1 +
n
n
φ
ψ
φ
φ
u≥0
+
−
i
i
1
2
i
i
−
−1
⎪⎪
2
=⎨
⎪φ n − 1 ψ − φ n − φ n
u<0
i+1
⎪⎩ i + 1 2 i + 3 2 i + 2
(
)
ここで、ψはリミッターで近傍のφの分布によって値を変化させる
n
n
−
φ
φ
+
i
+
1
i
ψ i+− 1 2 = ψ i+− 1 2 (ri+− 1 2 ) ただし ri − 1 2 = n
φi − φin− 1
ψ i−+ 3 2 = ψ i−+ 3 2 (ri−+ 3 2 )
φin+ 1 − φin
−
ri + 3 2 = n
φi + 2 − φin+ 1
TVD条件を満たす代表的なリミッターには、以下のminmod関数がある
ψ i+− 1 2 = min mod (1 , ri+− 1 2 )
⎧x
⎪
min mod ( x , y ) = ⎨ y
⎪0
⎩
x ≤ y and xy ≥ 0
x > y and xy ≥ 0
xy < 0
最新の空間差分スキーム
- TVD (Total Variation Diminishing) (2/2)φin+ 1 2
(
)
⎧ n 1 +
n
n
u≥0
φ
ψ
φ
φ
+
−
i −1
ψ i+− 1 2 = ψ i+− 1 2 ri+− 1 2
⎪⎪ i 2 i − 1 2 i
=⎨
ψ i−+ 3 2 = ψ i−+ 3 2 ri−+ 3 2
⎪φ n − 1 ψ − φ n − φ n
u<0
i+1
⎪⎩ i + 1 2 i + 3 2 i + 2
(
(
(
)
ψ i+− 1 2 = min mod (1 , ri+− 1 2 )
⎧x
⎪
min mod ( x , y ) = ⎨ y
⎪0
⎩
⎧1
⎪
minmod (1, r ) = ⎨r
⎪0
⎩
)
)
+
i −1 2
φin+ 1 − φin
= n
φi − φin− 1
−
i+3 2
φin+ 1 − φin
= n
φi + 2 − φin+ 1
r
r
x ≤ y and xy ≥ 0
x > y and xy ≥ 0
xy < 0
1≤r
0≤r <1
r <0
これにより、３種類の差分スキームを状況に応じて使い分けることができる。
・ 1<rの場合：変数の勾配が風下方向の増加 → ２次風上差分を適用
・ 0<r<1の場合：変数の勾配が風下方向に減少 → ２次中心差分を適用
・ r<0の場合：変数の勾配が逆転 → １次風上差分を適用
最新の空間差分スキーム
- CIP (Cubic-interpolated pseudo particle) 物理量の空間１階微分を変数として持ち、その移流方程式を同時に解く。
∂φ
∂φ
= −u
∂t
∂x
∂ ⎛ ∂φ ⎞
∂ ⎛ ∂φ ⎞
⎜
⎟ = −u ⎜
⎟
∂ t ⎝ ∂x ⎠
∂x ⎝ ∂x ⎠
さらに、φの空間分布を３次の補間式で近似する。
φ ( x ) = a0 + a1 x + a 2 x 2 + a 3 x 3
この補間式の４つの係数は、位置iと位置i-1での変数と空間微分値によって、
次式によって決定する
n
a0 = φ i
⎛ ∂φ ⎞
⎟
⎝ ∂x ⎠ i
φin+ 1 = φin − Δtu⎜
a1 = (φ , x )i
[2(φ , x ) + (φ , x ) ] − 3 [φ
=
i
i −1
i − φi − 1 ]
a2
Δx 2
Δx
[(φ , x )i + (φ , x )i −1 ] − 2 [φi − φi −1 ]
a3 =
Δx 2
Δx 3
n+ 1
⎛ ∂φ ⎞
⎜
⎟
⎝ ∂x ⎠ i
⎧ ⎡ ⎛ ∂φ ⎞ n ⎛ ∂φ ⎞ n ⎤
⎪ ⎢ 2⎜
⎟ +⎜
⎟ ⎥
n
∂
∂
x
x
⎝
⎠
⎝
⎠ i − 1 ⎥⎦
⎪
φin − φin− 1
∂
φ
⎢
i
⎞
⎛
⎣
=⎜
−3
⎟ − Δtu 2 ⎨
Δ
x
Δx 2
⎝ ∂x ⎠ i
⎪
⎪
⎩
[
⎫
⎪
⎪
⎬
⎪
⎪
⎭
]
熱伝導方程式(放物型)の差分スキーム
熱伝導方程式の陽解法(Explicit scheme)
熱伝導方程式
∂u
∂ 2u
=α 2
∂t
∂x
差分方程式
( )
uin + 1 − uin
uin+ 1 − 2uin + uin−1
2
(
)
=α
+
Δ
+
Δ
O
t
O
x
Δt
Δx 2
uin+1 について解くと、以下のFTCSスキームが得られる。
Δt
n+1
n
n
n
n
d =α
ui = ui + d ui + 1 − 2ui + ui −1
Δx 2
FTCSスキーム
d：拡散数（diffusion number）
(
)
熱伝導方程式の陽解法
Δx
差分方程式
(
uin + 1 = uin + d uin+ 1 − 2uin + uin−1
Δt
d =α
Δx 2
陽解法によるプログラム例
1590
1600
1610
1620
)
n+1
Δt
n
n-1
t
0
i-1
i i+1
*FTCS
FOR I=2 TO IM1
U(I)=UN(I)+D*(UN(I+1)-2*UN(I)+UN(I-1))
NEXT I:RETURN
x
熱伝導方程式の陰解法 (Implicit scheme)
熱伝導方程式：
∂u
∂ 2u
=α 2
∂t
∂x
時間に後退差分近似、空間に中心差分近似
した差分方程式：
( )
uin + 1 − uin
uin++11 − 2uin + 1 + uin−+11
2
(
)
O
t
O
x
=α
+
Δ
+
Δ
Δt
Δx 2
uin−+11 , uin+1 , uin++11 を変数値として差分スキームを求める。
n+1
i −1
− du
+ (1 + 2d )u
n+1
i
n+1
i +1
− du
=u
n
i
Δt
d =α
Δx 2
熱伝導方程式の陰解法
差分方程式：
− duin−+11 + (1 + 2d )uin+1 − duin++11 = uin
(例) i=1~5を考える。i=1, i=5を
境界条件を満足する境界点と
し、i=2, 3, 4を内点とする。
(1 + 2d )u2n+1 − du3n+1 = u2n + du1n+1
− du2n+1 + (1 + 2d )u3n +1 − du4n +1 = u3n
− du3n+1 + (1 + 2d )u4n+1 = u4n + du5n+1
N元連立一次方程式
境界
tn
n+1
n
i
1
2
3
4
5
熱伝導方程式の陰解法
三重対角行列
⎛ (1 + 2d )
⎜
⎜
⎜
⎝
−d
−d
(1 + 2d )
0
−d
n
n+1
n+1
⎛
⎞
⎛
⎞
+
u
du
u
0 ⎞⎜ 2 ⎟ ⎜ 2
1
⎟
⎟ n+1
⎟
− d ⎟⎜ u3 ⎟ = ⎜ u3n
(1 + 2d )⎟⎠⎜⎜⎝ u4n+1 ⎟⎟⎠ ⎜⎜⎝ u4n + du5n+1 ⎟⎟⎠
陰解法によるプログラム例(行列の設定)
2020
2030
2040
2050
2060
*SPMTRX
FOR I=1 TO IBAR
A(I,I)=1+2*BETA:NEXT I
FOR I=1 TO IBAR-1:A(I,I+1)=-BETA:NEXT I
FOR I=1 TO IBAR-1:A(I+1,I)=-BETA:NEXT I:RETURN
熱伝導方程式の陰解法
陰解法によるプログラム例（三重対角行列の解法）
1580
1590
1591
1600
1610
1620
1630
1640
1650
1660
1670
1680
1690
1700
1710
'*Gauss-Jordan eliminator
*GSSJOD
' *forward elimination
FOR K=1 TO IBAR
D1(K)=1/A(K,K)
FOR I=1 TO IBAR
IF I=K THEN GOTO 1700
M(I,K)=A(I,K)*D1(K)
FOR J=K+1 TO JBAR
A(I,J)=A(I,J)-M(I,K)*A(K,J)
NEXT J
B(I)=B(I)-M(I,K)*B(K)
A(I,K)=0
NEXT I
NEXT K
1720
1730
1740
1760
'*back substitution
FOR I=1 TO IBAR
U(I)=B(I)*D1(I)
NEXT I
歴史的な経緯 - 初期の差分法
Richardson(1881-1953) の方法

差分の考え方は１８６０年代、既にあった。しかしながら、
放物型方程式に対して初めて差分法を適用したのは、
Richardon(1881-1953)で、1910年のことである。
Richardon は１次元熱伝導方程式対して、時間と空間微分
に中心差分を用いることにより、以下の差分式を導いた。
∂u
∂ 2u
=α 2
∂t
∂x
uin + 1 − uin − 1
uin− 1 − 2 uin + uin+ 1
=α
2 Δt
Δx 2
n+ 1
i
u

n− 1
i
=u
(
+d u
n
i −1
− 2u + u
n
i
n
i+1
1
uin +1
d
-2d
uin−1
)
d =α
Δt
Δx 2
d
uin+1
uin
1
uin −1
残念ながら、このスキームは、拡散数dの値によらず不安定となるため、正しい解
が得られないスキームであった。
近年、このRichardonスキームは不安定なスキームの代表例として扱われることが
あるが、初めて実用的な計算に差分法を適用しようとした試みとして高く評価され
るものである。
Richardsonの夢

リチャードソンは気象学に興味を持ち、微分方程式を離散
化して数値的に解くという現在使われている気象予測の原
型ともいえる数値予報システムを提案した。
そのために、大戦中に集めたデータを用いた計算を実際に
行ったが、計算をすべて手で行ったためわずか6時間の予
報に2ヶ月かかり、しかも数値の処理に問題があったため、
予報は失敗に終わった。
1922年に著書の中で彼は、64,000人の計算者を巨大な
ホールに集めて指揮者の元で整然と計算を行えば実際の
天候の変化と同じくらいの速さで予報が行えると見積もっ
た。この構想は「リチャードソンの夢」と呼ばれた。
またリチャードソンは乱気流に関する実験も数多く行ってい
る。乱流に関する無次元数「リチャードソン数」は彼の名前
にちなんでいる。
デュフォール・フランケル法(DF法)

Richardson法は、拡散数に関わらず不安定な差分スキームであったが、
DufortとFrankelは、１９５３年、空間の２階微分に対する差分表示を時間平
均値とする差分スキームを導いた。
このスキームは、拡散方程式に対するLeapfrog法(あるいはPyramid法)とも
呼ばれているが、対流方程式に対するLeapfrog法と異なり、安定なスキー
ムとなっている。（→ 安定解析により証明）
uin− 1 − 2 uin + uin+ 1
∂ 2u
=α
∂x 2
Δx 2
uin +1
uin + 1 + uin − 1
u =
2
n
i
uin+1
uin−1
(
)
uin+ 1 − uin − 1
uin− 1 − uin+ 1 + uin− 1 + uin+ 1
=α
2
2 Δt
Δx 2
uin+ 1 =
[ (
)
1
2 d uin− 1 + uin+ 1 + (1 − 2 d )uin− 1
1 + 2d
uin −1
]
d =α
Δt
Δx 2
アダムス・バシュフォース法(AB法)
⎡ uin− 1 − 2 uin + uin+ 1
uin + 1 − uin
uin−−11 − 2 uin− 1 + uin+−11 ⎤
= α ⎢β
+ (1− β )
⎥
2
2
Δt
Δ
x
Δ
x
⎣
⎦
3
β=
2
3
1
n+ 1
n
n
n
ui = dui − 1 + (1 − 2 d )ui + dui + 1 − duin−−11 + (1 − 2 d )uin − 1 + duin+−11
2
2
{
} {
}
n+1
Δt
d =α
Δx 2
n
n-1
i-1
i
i+1
クランク・ニコルソン法(CN法) ← 陰解法

CN法は、空間微分に対してnとn-1の２つの時間レベルの算術平均を
取っている。
得られる漸化式は、各種の行列演算を用いて解かれる。
∂u
∂ 2u
=α 2
∂t
∂x
f =
n
n
i −1
u
uin − uin− 1 1 n
=
f + f n− 1
2
Δt
{
n
− 2u + u
Δx
n
i
2
n
i+1
}
n-1
i-1
i
i+1
− duin− 1 + 2(1 + d )uin− 1 − duin+ 1 = duin−−11 + 2(1 − d )uin− 1 + duin+−11 d = α
Δt
Δx 2
Poisson方程式(楕円型)の差分スキーム
（1/6）中心差分スキームの適用
Poisson方程式
∂ 2u ∂ 2u
+ 2 = −g
2
∂y
∂x
空間に中心差分近似を用いた差分方程式
ui +1, j − 2ui , j + ui −1, j
Δx
2
+
ui , j +1 − 2ui , j + ui , j −1
Δy
2
= − gi , j
Δx=Δyとする。
− ui +1, j − ui −1, j + 4ui , j − ui , j +1 − ui , j −1 = Δx 2 g i , j
（2/6） Natural ordering
Natural ordering:
u1,1, u1,2,……, u3,3をu1, u2, ……, u9と番号を付け直す
y
y
21
22
23
24
7
8
9
4
5
6
1
2
3
11
12
13
25
(4,0) (4,1) (4,2) (4,3) (4,4)
17
(3,0) (3,1) (3,2) (3,3) (3,4)
16
(2,0) (2,1) (2,2) (2,3) (2,4)
15
(1,0) (1,1) (1,2) (1,3) (1,4)
(0,0) (0,1) (0,2) (0,3) (0,4)
x
10
20
19
18
14
x
（3/6）行列式
Natural ordering によって差分式は以下のように
表すことができる
−1
⎛ 4 −1
⎞⎛ u1 ⎞ ⎛ Δx 2 g1 + u11 + u15 ⎞ ⎛ b1 ⎞
⎟ ⎜ ⎟
⎜
⎟⎜ ⎟ ⎜ 2
u
−1
⎟ ⎜ b2 ⎟
⎜−1 4 −1
⎟⎜ 2 ⎟ ⎜ Δx g 2 + u12
⎜
⎟⎜ u ⎟ ⎜ Δ x 2 g + u + u ⎟ ⎜ b ⎟
−1 4
−1
3
13
18 ⎟
⎜ 3⎟
⎜
⎟⎜ 3 ⎟ ⎜
⎟ ⎜ b4 ⎟
4 −1
−1
⎜−1
⎟⎜ u4 ⎟ ⎜ Δx 2 g 4 + u16
⎟ ⎜ ⎟
⎜
⎟⎜ u ⎟ = ⎜ 2
1
1
4
1
1
−
−
−
−
x
g
Δ
5
⎟ ≡ ⎜ b5 ⎟
⎜
⎟⎜ 5 ⎟ ⎜
⎟ ⎜ b6 ⎟
⎜
−1
−1 4
− 1 ⎟⎜ u6 ⎟ ⎜ Δx 2 g6 + u19
⎟ ⎜ ⎟
⎜
⎟⎜ ⎟ ⎜ 2
4 −1
−1
⎜
⎟⎜ u7 ⎟ ⎜ Δx g7 + u17 + u22 ⎟ ⎜ b7 ⎟
⎟ ⎜b ⎟
⎜
−1
− 1 4 − 1 ⎟⎜ u8 ⎟ ⎜ Δx 2 g8 + u23
⎟ ⎜⎜ 8 ⎟⎟
⎟⎟⎜⎜ ⎟⎟ ⎜⎜ 2
⎜⎜
−1
− 1 4 ⎠⎝ u9 ⎠ ⎝ Δx g9 + u20 + u24 ⎟⎠ ⎝ b9 ⎠
⎝
（4/6）行列式の要素～三重対角行列
U1 = (u1
u2
B1 = (b1 b2
u3 ) T , U 2 = (u4
b3 ) T ,
B2 = (b4
u5
b5
u6 ) T , U 3 = (u7
b6 ) T ,
B3 = (b7
u8
b8
u9 ) T
b9 ) T
以上の置き換えによって、行列式は、以下の形に書き換わる
⎛ A11
⎜
⎜ A21
⎜
⎝
ここで
A12
A22
A32
⎞⎛ U 1 ⎞ ⎛ B1 ⎞
⎟⎜ ⎟ ⎜ ⎟
A23 ⎟⎜ U 2 ⎟ = ⎜ B2 ⎟
A33 ⎟⎠⎜⎝ U 3 ⎟⎠ ⎜⎝ B3 ⎟⎠
⎛ 4 −1
⎞
⎜
⎟
1
4
1
Aii = ⎜ −
− ⎟
⎜
⎟
1
4
−
⎝
⎠
⎛−1
⎞
⎜
⎟
−1
Aij = A ji = ⎜
⎟
⎜
⎟
1
−
⎝
⎠
(i = 1,2,3)
(i , j = 1,2,3 : i ≠ j )
（5/6）ブロック三重対角行列の解法
係数行列の設定：
A行列の設定：
1710
1720
1740
1750
1760
1770
1780
1790
1800
1810
1820
'*specification of A
1830 FOR J0=1 TO JBAR-1
*MTRX
1840 FOR I=1 TO IBAR:FOR J=1 TO JBAR
FOR I=1 TO IBAR:AII(I,I)=4:NEXT I
1850 IP=I+JBAR*J0: JP=J+JBAR*(J0-1)
FOR I=2 TO IBAR:AII(I,I-1)=-1:NEXT I
1860 A(IP,JP)=AIJ(I,J):NEXT J:NEXT I
FOR I=1 TO IBAR-1:AII(I,I+1)=-1:NEXT I 1870 NEXT J0
FOR I=1 TO IBAR:AIJ(I,I)=-1:NEXT I
1880 FOR J0=1 TO JBAR-1
FOR J0=1 TO JBAR
1890 FOR I=1 TO IBAR:FOR J=1 TO JBAR
FOR I=1 TO IBAR:FOR J=1 TO JBAR
1900 IP=I+JBAR*(J0-1):JP=J+JBAR*J0
IP=I+JBAR*(J0-1):JP=J+JBAR*(J0-1)
1910 A(IP,JP)=AIJ(I,J):NEXT J:NEXT I
A(IP,JP)=AII(I,J):NEXT J:NEXT I
1920 NEXT J0:RETURN
NEXT J0
b行列の設定：
1940
1950
1960
1970
'*specification of b
*RHS
FOR I=1 TO IJBAR:B(I)=G0*DXDY:NEXT I:RETURN
'FOR I=IBAR TO IJBAR STEP IBAR:B(I)=U0:NEXT I:RETURN
（6/6）ブロック三重対角行列の解法
ガウスの消去法による行列の解法
1440
1450
1470
1480
1490
1500
1510
1520
1530
1540
1550
1560
1580
1590
'*Gaussian eliminator for A*u=b
*GAUSS
FOR K=1 TO IJBAR
D1(K)=1/A(K,K)
FOR I=K+1 TO IJBAR
M(I,K)=A(I,K)*D1(K)
FOR J=K+1 TO IJBAR
A(I,J)=A(I,J)-M(I,K)*A(K,J)
NEXT J
B(I)=B(I)-M(I,K)*B(K)
A(I,K)=0
NEXT I
NEXT K
D1(IJBAR)=1/A(IJBAR,IJBAR)
1600
1610
1620
1630
1640
1650
1660
1690
'*backward substitution
U(IJBAR)=B(IJBAR)*D1(IJBAR)
FOR I= 2 TO IJBAR
IR1=IJBAR-I+1
SUM=0:FOR J=IR1+1 TO IJBAR
SUM=SUM+A(IR1,J)*U(J):NEXT J
U(IR1)=(B(IR1)-SUM)/A(IR1,IR1)
NEXT I
5点差分スキーム
(Five-term recurrence scheme)
Δx=Δyとした場合のPoisson方程式の差分近似式
− ui +1, j − ui −1, j + 4ui , j − ui , j +1 − ui , j −1 = Δx 2 g i , j
ui,jについて整理すると、
ui , j =
(
1
ui + 1, j + ui −1, j + ui , j + 1 + ui , j −1 + Δx 2 gi , j
4
)
変数、ui+1,j, ui-1,j, ui,j+1, ui,j-1 が既知の場合、ui,j の値が
求まることを示している。
表計算ソフトによる熱伝導解析(1/2)
100
51.07237
31.6532
18.33154
9.599884
4.507474
1.896977
0.716155
0.227735
0
100
64.4262
39.86327
22.32332
11.13954
4.936661
1.953679
0.697562
0.213375
0
100
66.64282
39.9896
21.16744
9.804685
3.986152
1.448326
0.483257
0.142508
0
100
65.19423
36.95083
17.95006
7.44554
2.683597
0.878807
0.27722
0.079438
0
100
61.99906
32.13504
13.67962
4.818129
1.479724
0.45038
0.136799
0.038443
0
100
57.19167
25.61034
8.787319
2.370757
0.668902
0.195639
0.058171
0.016165
0
100
80
80-100
60-80
40-60
20-40
0-20
60
40
20
9
7
5
S5
3
0
S9
1
100
0
0
0
0
0
0
0
0
0
S1
100
49.72089
17.04673
3.538683
0.884245
0.241968
0.069671
0.020536
0.00568
0
100
37.43022
4.261682
0.884671
0.221061
0.060492
0.017418
0.005134
0.00142
0
100
0
0
0
0
0
0
0
0
0
表計算ソフトによる熱伝導解析(2/2)
20
20
20
20
20
20
20
20
20
20
20
20.36674
20.76934
21.24008
21.79738
22.41812
22.97566
23.13608
22.29624
20
20
20.69783
21.47087
22.39394
23.53172
24.89976
26.34872
27.27259
26.04895
20
20
20.954
22.02284
23.33364
25.03634
27.30095
30.24724
33.55684
34.62707
20
20
20
21.09569 21.09575
22.33339 22.3335
23.88207 23.88221
25.97965 25.9798
29.02103 29.02117
33.78289 33.78301
42.08071 42.08079
58.90258 58.90262
100
100
加熱部
20
20.95416
22.02314
23.33402
25.03674
27.30133
30.24755
33.55705
34.62718
20
20
20.69803
21.47124
22.39441
23.53221
24.90022
26.3491
27.27286
26.04909
20
100
90
80
70
60
50
40
30
20
S9
S7
S5
10
0
1
2
3
4
5
S3
6
7
8
9
10
S1
90-100
80-90
70-80
60-70
50-60
40-50
30-40
20-30
10-20
0-10
20
20.3669
20.76962
21.24042
21.79775
22.41846
22.97595
23.13628
22.29634
20
20
20
20
20
20
20
20
20
20
20
差分法の特徴
(1/5)微係数を差分商に置き換える
対流方程式の偏微分方程式
∂u
∂u
+c
=0
∂t
∂x
差分方程式
uin+1 − uin
uin − uin−1
+c
=0
Δt
Δx
微係数を差分商に置き換える操作は数学的表現
として、元の方程式の形が大きく崩れないで、差分
方程式に引き継がれる。
（2/5）物理の過程を忠実に再現
対流方程式の差分方程式：
(
)
uin+1 = uin − C uin − uin−1 , C = c
Δt
Δx
初期値問題において、時間の進展を一つ一つのΔtの
増分として現象を追跡している。
Laplace方程式の差分式 (gi,j=0とし、u：物質の濃度）
ui , j =
(
1
ui + 1, j + ui −1, j + ui , j + 1 + ui , j −1
4
)
周囲の4点の平均として、注目している点の濃度が決定され
る。という物理学における拡散の過程までもが理解される。
（3/5）数学的な意味付けが容易
偏微分方程式の従属変数をTaylor展開して差分方程式に変
換したとみなすと、
打切り誤差、境界条件の誤差の評価が容易となり、高次項ま
で考慮した差分法の数学的な意味付けも容易に得られる。
（4/5）汎関数の存在を必要としない
スカラー関数u(x,y)に対するDiirichlet積分は、対象領域をD
として
⎡ ⎛ ∂u ⎞ 2 ⎛ ∂u ⎞ 2 ⎤
I = ∫D ⎢⎜ ⎟ + ⎜⎜ ⎟⎟ ⎥dxdy
⎢⎣⎝ ∂x ⎠ ⎝ ∂y ⎠ ⎥⎦
で定義される。この汎関数を最小とする問題のEularLagrangeの方程式は、次のLaplace方程式である。
∂ 2u
∂ 2u
+c 2 =0
2
∂y
∂x
差分法では、対象としている偏微分方程式についての汎関
数が存在するか否かを気にかける必要がない。
（5/5）検証が容易
差分スキームによる計算アルゴリズムは、物理的直感が働
きやすく、定式化と変換の次にくる計算過程の設計や検証が
容易である。例えば、
(
)
uin+1 = uin − C uin − uin−1 , C = c
Δt
Δx
においてC=1としたとき、次の厳密な差分解
uin+1 = uin−1
が得られる。この場合にはΔtが1つ増加すると、点(i-1,n)の変
数uがxの正方向にΔxだけ移動して点(i,n+1)に進行するとい
う有向性の意味が明確になる。
課題提出予定(案)
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
( 9/ 2) 数値シミュレーションの手続き (テキスト第１章)
( 9/ 9) 課題(1)提出
( 9/16) 休講
( 9/30) 課題(2)提出
(10/ 7) 課題(3)提出
(10/21) 課題(4)提出
(10/28) 課題(5)提出
(11/ 4) 課題(6)提出
(11/11) 課題(7)提出
(11/18) (課題by金子暁子先生)の提出
(11/25) (数値解析by渡辺正先生)の提出
提出課題(1)

各自が行ったことのある数値解析（流体解析に限定
しない）について、その対象となった物理、基礎方程
式、境界条件、数値解析手法、解析結果、（もしあ
れば）実験との比較、等に関して、簡潔にレポートに
纏めなさい。
もし、数値解析の経験がない場合には、数値解析
に関して考えるところをレポートに纏めなさい。
提出期限は、次回（9/9）までとする。
提出課題(2)

熱伝導方程式
∂ 2ϕ
∂ϕ
=α 2
∂t
∂x

(α > 0)
に、次の関数を代入すると、
⎛ 1
⎞
∫ udx ⎟
⎠
⎝ 2α
ϕ = exp⎜ −

以下の式が導かれることを示せ。
∂ 2u
∂u
∂u
=α 2
+u
∂x
∂x
∂t

提出期限は、次々回（9/30）までとする。
提出課題(3)

以下のテンソル表現で記載されているNavierStokes方程式の各項の成分を書き下しなさい。
r
ρ

r
Dv
= −∇p − ∇ ⋅ π + ρg
Dt
ただし、 vr = (v x , v y , v z )
r
r
r
r r ∂v r r
Dv ∂v
=
+ ∇ ⋅ (v v ) =
+ v ⋅ ∇v
Dt ∂t
∂t
[
r
rT
∇ ⋅ π = −μ∇ ⋅ ∇v + (∇v )

]
提出期限は、次次々回（10/7）までとする。
対流方程式の数値解の挙動
～スキームの選択ーモデル問題３－

以下の対流方程式を
∂u
∂u
+c
=0
∂t
∂x

以下の境界条件と初期条件の下で、
⎧0 ( 0 ≤ x ≤ 1 / 5)
⎪
u( x ,0) = F ( x ) = ⎨1 (1 / 5 ≤ x ≤ 2 / 5)
⎪0 ( 2 / 5 ≤ x ≤ 1)
⎩

u(0, t ) = g0
g0 = 0
∂u(1, t )
=0
∂x
以下の各スキームを用いて数値的に解析しなさい。
•
•
•
FTCS
蛙とび法
風上差分
Program: ADVEC
提出課題(4)

前ページに指定する微分方程式を、与えら
れた境界条件の下で、指定する差分スキー
ムを用いて数値的に解析しなさい。
提出期日： 10月21日（木）
提出物は、
・内容の説明したレポート
・プログラムリスト
・解析結果（数値出力、作図結果）
以上の全てを、A4レポート用紙ならびに電
子媒体(CD)の両方により提出すること。
提出課題(5)
以下のPoisson方程式を、中心差分近似を用いた差分方
程式とした上で、
∂ 2u ∂ 2u
+ 2 = −g
2
∂y
∂x
５×５に分割された空間に対して、周囲境界を既知のもの
とするとき得られる９元連立一次方程式を、Natural
orderingによって書き下しなさい。
提出期日：10月28日（木）

Download Report