偏微分方程式の数値解法

地球惑星科学実習 B-3: 偏微分方程式の数値解法
岩山 隆寛 ∗†
2015 年 7 月 3, 10, 17, 24 日
1 はじめに
地球惑星科学の諸現象は, 偏微分方程式の形で書かれることが多い. 線形の偏微分方程
式であれば多くの場合, 解析的に（手計算で）解くことができるが, 非線形の偏微分方程式
や, 境界条件が複雑な場合には, 解析的に解くことは困難である. このようなときには計算
機を用いて, その偏微分方程式を解く. ここでは, 偏微分方程式の数値解法の入門として,
解析的に解くことができるが, 拡散方程式や波動方程式などの簡単な偏微分方程式を計算
機を用いて解いてみる.
本実習では, 簡単化のために空間 1 次元とし, 以下の偏微分方程式を取り扱う.
1. 拡散方程式:
∂u
∂2u
= ν 2,
∂t
∂x
(1)
∂u
∂u
+c
= 0,
∂t
∂x
(2)
ここで, ν は正の定数である.
2. 線形移流方程式:
ここで, c は定数である.
∗
†
神戸大学 大学院理学研究科 惑星学専攻. e-mail: [email protected]
本実習の資料の作成に際して, 名古屋工業大学 大学院工学研究科 渡邊威 准教授と神戸大学 大学院理学
研究科研究員 村上真也 博士の協力を頂きました.
1
3. 線形移流拡散方程式:
∂u
∂u
∂2u
+c
= ν 2,
∂t
∂x
∂x
(3)
ここで, c は定数, ν は正の定数である.
4. Burgers 方程式:
∂u
∂2u
∂u
+u
= ν 2,
∂t
∂x
∂x
(4)
ここで, ν は正の定数である.
これらの方程式は, 流体力学の基礎方程式である (1 次元の）Navier–Stokes 方程式,
∂u
1 ∂p
∂2u
∂u
+u
=−
+ ν 2,
∂t
∂x
ρ ∂x
∂x
(5)
を想定している. ここで, ν は正の定数である.
2 gnuplot による png ファイルの作成と, アニメーションの
作り方
偏微分方程式の解は, 物理量の空間分布が時間発展するので, 解をアニメーションで表
示すると見やすいであろう. ここでは, 与えられた関数を gnuplot を使用して png ファイ
ルを作り, そのファイルから（パラパラ漫画式の）アニメーションを作成する方法を解説
する. ここで解説する方法は, 計算結果からアニメーションを作成する方法にも応用可能
である.
2.1 スクリプトファイルを利用した png ファイルの作成方法
gnuplot を利用して多くの図を作成するときに, gnuplot のコマンドをいちいち打つこ
とは効率が悪い. gnuplot のコマンドをファイル（スクリプトファイルと呼ばれる）に記
述しておき, それを gnuplot で読み込んで図を作成する方法を解説する.
1. test.plt という名前のスクリプトファイルを作成する.
2. 次の内容を, 適当なエディタを使用して, test.plt に記述する.
2
set xrange[0:2*pi]
set yrange[-1:1]
set xlabel ’x’
set ylabel ’u(x,t)’
plot sin(x) title ’sin(x)’
set term png
set output ’sincurve.png’
replot
3. gnuplot を起動して, 以下のコマンドを打ち, sin(x) が画面に表示され, さらに
sincurve.png ファイルができ, そのファイルに表示された図と同じ図 (図 1 参照)
が書き込まれていることを確認する.
gnuplot> cd ’test.plt を保存したディレクト名’
gnuplot> load ’test.plt’
1
sin(x)
u(x,t)
0.5
0
-0.5
-1
0
1
2
3
x
4
図 1 sin(x) のグラフ.
3
5
6
2.2 gif アニメーションの作成方法：その１
1. 以下の内容のスクリプトファイル（test2.plt) を作成し, gunplot で実行しなさい.*1
set xrange[0:2*pi]
set yrange[-1:1]
set xlabel ’x’
set ylabel ’u(x,t)’
set term png
set size 1, 1
set output ’sincurve_000.png’
plot sin(x) title ’t=0’
set output ’sincurve_001.png’
plot sin(x+pi/8.) title ’t=1’
set output ’sincurve_002.png’
plot sin(x+2.*pi/8.) title ’t=2’
.
.
.
set output ’sincurve_015.png’
plot sin(x+15.*pi/8.) title ’t=15’
2. コマンドプロンプトから ImageMagic の convert コマンドを使用して, png ファイ
ルを結合して gif アニメーションを作成する.
> convert sincurve_???.png sincurve.gif
3. sincurve.gif を適当なブラウザで開くと, ファイル名順に結合された gif アニメー
ションを見ることができる.
*1
eps ファイルでアニメーションを作ると, 前の時刻の画像が透けて見えるため, png ファイルでアニメー
ションを作ることにした.
4
2.3 gif アニメーションの作成方法：その２
fortran によって計算結果をファイルに書き出し, そのファイルに書かれているデータ
を gnuplot を用いて作図し, png ファイルとして保存する. さらに, convert コマンドを用
いて, png ファイルを結合してアニメーションを作成する方法の例を示す.
1. 次の fotran プログラム (test 1.f90) を解読しなさい.
! サンプルプログラム
! test_1.f90
! 作成者 渡邊威 (名古屋工業大学 工学研究科 シミュレーション工学専攻)
! 改変: 2015.06.23 岩山隆寛
!--------------------------------------------------------------------implicit none
integer
n, i, j, it, nsteps, nout
real*8
t, h, tout, pi, nu, Lx, l, a
parameter (pi=3.141592653589d0, Lx=2.0d0*pi, nu=2.0d-3)
parameter (n=256, t=40.0d0, h=1.0d-2, tout=1.0d+0)
real*8
u(0:n+1), st(n,0:int(t/tout))
l=Lx/dble(n)
nsteps=int(t/h)
nout=int(tout/h)
it=0
! initial condition
do i=1,n
u(i)=dsin(dble(i)*l)
enddo
! 初期条件の保存
do i=1,n
st(i,it)=u(i)
enddo
5
! 時間発展
do j=1,nsteps
do i=1,n
u(i)=dsin(dble(i)*l+dble(j)*pi/dble(nsteps))
enddo
! データの保存
if (mod(j,nout)==0) then
it=it+1
do i=1,n
st(i,it)=u(i)
enddo
endif
enddo
! ファイルへの書き出し
do it=0, int(t/tout)
do i=1,n
write(it+100,’(2f10.3)’) i*l, st(i,it)
enddo
enddo
stop
end
2. test 1.f90 を実行し, fort.100∼fort.140 の 41 個のファイルができていることを確
認しなさい.
> gfortran test_1.f90
> ./a.out
3. gnuplot のスクリプトファイル (test 1.plt)*2 を実行して, fort.100∼fort.140 の
ファイルの中身を作図しなさい. (mov100.png∼mov140.png ができる.)
*2
これは以下のディレクトリにあるので, 各自自分の作業しているディレクトリにコピーして使用してくだ
さい.
6
gnuplot> load ’test_1.plt’
4. convert コマンドを打って, mov100.png∼mov140.png を結合して一つのファイル
move.gif を作成する.
> convert mov???.png move.gif
2.4 課題
test 1.f90 を自分で改良して, diffusion thoery.f とし, 無限領域内の 1 次元拡散方程式
の基本解,
{
x2
√
exp −
4ν(t + t0 )
4πν(t + t0 )
1
}
(6)
の発展をアニメーションで示しなさい. 拡散係数 ν や時間間隔, 等は自分で適当な値を設
定しなさい.
作 図 し た 関 数 に 関 す る 説 明 と gif フ ァ イ ル を 添 付 フ ァ イ ル と し て e-mail で
[email protected] 宛に送信してください.
3 拡散方程式の数値解法
1 次元拡散方程式,
∂u
∂2u
= ν 2,
∂t
∂x
(7)
を数値的に解く. ここで, ν は正の定数である.
3.1 解説
• 時間と空間の離散化:
拡散方程式 (7) を解く領域は, 0 ≤ x ≤ L とする. この区間を N 等分する. そこ
で, 座標変数 x を正の整数 i を用いて,
x = xi = i∆x,
\home3\iwayama\exp_15\pde\exp_1\test_1.plt
7
∆x =
L
,
N
(8)
と表せる. また, 時間も数値計算では離散化して扱うので, j を正の整数として,
t = tj = j∆t,
(9)
と表す. (7) の解は, x, t に依存する量なので, したがって, 数値計算では (7) の解
は, i, j に依存することになる:
u(x, y) → u(xi , tj )
(10)
• 空間微分の表現の仕方:
(7) は時間と空間に関する微分を含むが, 数値計算はでは微分は近似的にしか取り
扱えない. ここでは空間微分の近似の仕方について述べる.
xi における ∂ 2 u/∂x2 を xi−1 , xi , xi+1 の 3 点における u の値を使って表現
してみる. ∆x が微小のとき, u(xi±1 , tj ) を Taylor 展開すると
u(xi±1 , tj ) = u(xi , tj ) ±
∂u
1 ∂2u
∆x +
(∆x)2 + O(∆x3 ),
∂x
2 ∂x2
(11)
となる. ここで, O(∆x)3 は (∆x)3 の大きさの項を表す. したがって, 上式を
∂ 2 u/∂x2 について解いて
∂2u
u(xi+1 , tj ) + u(xi−1 , tj ) − 2u(xi , tj )
=
+ O(∆x2 )
2
∂x
(∆x)2
(12)
を得る.*3 以上から, (∆x)2 の精度で
∂2u
u(xi+1 , tj ) + u(xi−1 , tj ) − 2u(xi , tj )
=
2
∂x
(∆x)2
(13)
と近似する.
• 時間微分の表現の仕方:
時間微分に関しては, 常微分方程式の解法の際に採用した Euler 法を使用する. こ
のとき, ∆t の 1 次の精度で,
u(xi , tj+1 ) − u(xi , tj )
∂u
=
∂t
∆t
と近似する.
*3
O(∆x) 項は相殺されて, 誤差は O(∆x)2 程度の大きさになることに注意.
8
(14)
• まとめ:
以上の議論から, 拡散方程式 (7) を数値的に解く場合, 空間微分に関しては 2 次精
度, 時間微分に関しては 1 次精度で
u(xi , tj+1 ) − u(xi , tj )
u(xi+1 , tj ) + u(xi−1 , tj ) − 2u(xi , tj )
=ν
,
∆t
(∆x)2
(15)
もしくは,
u(xi , tj+1 ) = u(xi , tj ) + ν
∆t
{u(xi+1 , tj ) + u(xi−1 , tj ) − 2u(xi , tj )} ,(16)
(∆x)2
によって, 時刻 tj における u の分布から, tj+1 における u の分布が得られる.
3.2 課題
(16) に従って, 拡散方程式を数値的に解きなさい. ここで, 境界条件と初期条件は以下
のとおりとする.
境界条件：
周期境界条件,
u(x0 , tj ) = u(xN , tj ),
(17)
とする.
初期条件：
領域の中心を分布の中心とした Gauss 分布,
[
(
)2 ]
L
u(xi , t0 ) = exp −a x −
2
(18)
とする.
• 拡 散 方 程 式 を 数 値 シ ミ ュ レ ー シ ョ ン す る た め の サ ン プ ル プ ロ グ ラ ム (diffusion sample.f) を完成させる. 時間発展の様子をアニメーションにより観測す
る. gnuplot のスクリプトファイルは, test diffusion.plt を使用しなさい.
• 結果が妥当であるかどうかを検証するために, シミュレーションで得られた結果と
解析解を比較する. 解析解を計算するプログラムは, diffusion theory.f に用意され
ている.
• パラメータ ν を色々な値に変えてみてシミュレーションを行う. 結果をアニメー
ションにより観察し, 結果の意味をよく吟味する.
9
• 拡散方程式を数値的に解いた結果を図示した gif ファイルを提出しなさい.
– 計算条件（初期条件, 境界条件, 領域の大きさ, 解像度, ∆t) を明記する.
– 考察も付け加えること.(ν の値をかえたら. 超時間積分したら. など)
3.3 von Neumann の安定性解析
空間の差分の間隔 ∆t, 時間の差分間隔 ∆t が十分小さければ, 微分の差分による表現は
より正確になっていくであろう. しかしながら, 拡散係数 ν と ∆x, ∆t との間に以下の不
等式を満足する範囲内でしか, 拡散方程式は数値的に安定*4 には解けない:
(∆x)2
.
2ν
∆t ≤
(19)
(19) は von Neumann の安定性条件と呼ばれる. (19) によれば, ∆t, ν を固定したもとで
∆x を小さくしていくと数値計算は安定に実行できないことが分かる.
4 線形移流方程式の数値解法
1 次元線形移流方程式,
∂u
∂u
= −c ,
∂t
∂x
(20)
を数値的に解く. ここで, c は定数である.
4.1 移流方程式の解の性質
移流方程式 (20) は x − ct を独立変数とする任意の関数 f (x − ct) が解になっていると
いう性質がある:
u(x, t) = f (x − ct).
(21)
(21) は c が正のとき x の正の方向に速度 c で進行する解, c が負のとき x の負の方向に
速度 c で進行する解を表す.
*4
差分による近似によって生じた誤差が成長しないで解ける, という意味.
10
4.2 解説
• 時間と空間の離散化:
移流方程式 (20) を解く領域は, 0 ≤ x ≤ L とする. この区間を N 等分する. そこ
で, 座標変数 x を正の整数 i を用いて,
x = xi = i∆x,
∆x =
L
,
N
(22)
と表せる. また, 時間も数値計算では離散化して扱うので, j を正の整数として,
t = tj = j∆t,
(23)
と表す. (20) の解は, x, t に依存する量なので, したがって, 数値計算では (20) の
解は, i, j に依存することになる:
u(x, y) → u(xi , tj )
(24)
• 空間微分の表現の仕方:
(7) は時間と空間に関する微分を含むが, 数値計算はでは微分は近似的にしか取り
扱えない. ここでは空間微分の近似の仕方について述べる.
xi における ∂ 2 u/∂x2 を xi−1 , xi , xi+1 の 3 点における u の値を使って表現
してみる. ∆x が微小のとき, u(xi±1 , tj ) を Taylor 展開すると
u(xi±1 , tj ) = u(xi , tj ) ±
1 ∂2u
∂u
∆x +
(∆x)2 + O(∆x3 ),
2
∂x
2 ∂x
(25)
となる. ここで, O(∆x)3 は (∆x)3 の大きさの項を表す. したがって, 上式を
∂u/∂x について解くと, 3 つの表現が得られる:
1. 前進差分
∂u
u(xi+1 , tj ) − u(xi , tj )
=
+ O(∆x). ∂x
∆x
(26)
∂u
u(xi , tj ) − u(xi−1 , tj )
=
+ O(∆x). ∂x
∆x
(27)
2. 後進差分
11
3. 中央差分
∂u
u(xi+1 , tj ) − u(xi−1 , tj )
=
+ O(∆x2 ). ∂x
2∆x
(28)
• 時間微分の表現の仕方:
時間微分に関しては, 常微分方程式の解法の際に採用した Euler 法を使用する. こ
のとき, ∆t の 1 次の精度で,
∂u
u(xi , tj+1 ) − u(xi , tj )
=
∂t
∆t
(29)
と近似する.
• まとめ:
以上の議論から, 移流方程式 (20) を数値的に解く場合, 空間微分に関しては 1 次精
度, 時間微分に関しては 1 次精度で
1. 前進差分
u(xi , tj+1 ) = u(xi , tj ) +
c∆t
{u(xi+1 , tj ) − u(xi , tj )} . ∆x
(30)
u(xi , tj+1 ) = u(xi , tj ) +
c∆t
{u(xi , tj ) − u(xi−1 , tj )} . ∆x
(31)
2. 後進差分
空間微分に関しては 2 次精度, 時間微分に関しては 1 次精度で
3. 中央差分
u(xi , tj+1 ) = u(xi , tj ) +
c∆t
{u(xi+1 , tj ) − u(xi−1 , tj )} . 2∆x
(32)
によって, 時刻 tj における u の分布から, tj+1 における u の分布が得られる.
4.3 課題
(30)∼(32) に従って, 移流方程式を数値的に解きなさい. ここで, 境界条件と初期条件
は以下のとおりとする.
12
境界条件：
周期境界条件,
u(x0 , tj ) = u(xN , tj ),
(33)
とする.
初期条件：
波長 L/2, 振幅 1 の正弦波,
[
]
4π(x − ct)
u(xi , t0 ) = sin
,
L
(34)
とする.
その他の条件：
c = 0.1, L = 2π, N = 256, dt = 1.0 × 10−2 .
• 移流方程式を数値シミュレーションするためのプログラム (advection.f) を前回拡
散方程式を解くために使用したプログラム (diffusion sample.f) を書き換えて完成
させる. 時間発展の様子をアニメーションにより観測する. gnuplot のスクリプト
ファイルは, test diffusion.plt を適宜変更して使用しなさい.
• 結果が妥当であるかどうかを検証するために, シミュレーションで得られた結果と
解析解を比較する. 解析解を計算するプログラムは, diffusion theory.f に用意され
ている.
• パラメータ c の符号を変えてみてシミュレーションを行う. 結果をアニメーション
により観察し, これらの結果よりどの計算法が最も精度よく計算が行えているか吟
味する.
• 移流方程式を数値的に解いた結果を図示した gif ファイルを提出しなさい.
– 計算条件（初期条件, 境界条件, 領域の大きさ, 解像度, ∆t) を明記する.（注意:
プログラム上での表現でなく, 数学的な表現で書くこと.）
– 考察も付け加えること. (次の CFL 条件についても考慮して考察するとなおよ
い.)
• Courant-Friedrichs-Lewy(CFL) の条件について調べなさい.
4.4 Courant-Friedrichs-Lewy の条件
参考文献
• 川上 一郎, 「数値計算」, 第 6 章, 岩波書店, 1989 年.
13