実習の手引き
東京大学理学部・理学系研究科講義
地球物理数値解析
「差分法の基礎」
2014 年 4 月 27 日
横山央明
1
レポートについて
〆切： 実習実施日から 2 週間後
提出場所： 以下のどちらか。
(1) デジタル版の場合、横山 央明あて（[email protected]）、電子メールで。PDF ファイル
で送ること。
(2) 紙の場合。横山 央明居室（理学部 1 号館 805 号室）のドアポスト。自分用のコピーを控えとして保存
しておくこと。
注意：
(1) 形式は、A4 横書き で。表紙に以下を書くこと。
レポート課題名・問題出題のあった講義の日付・氏名・所属学部学科・学籍番号・提出日
(2) グラフのハードコピーを束ねただけのものはダメ。計算の設定や結果に対する考察を、一言でよいので
付けること。
返却： 地球惑星専攻宇宙惑星秘書室（理学部 1 号館 803 号室）にて、6 月 1 日から 6 月 30 日まで。それ以
後は廃棄する。
サンプル=プログラムのパッケージについて
この実習で用いるサンプル=プログラムのパッケージは、以下のウェブページからダウンロードすること。
gnuplot と Fortran コンパイラとを必要とする。自分の身のまわりにのこれらのソフトを使える環境がない
場合は、情報基盤センターの教育用計算機システムを申し込むとよい。
http://www-space.eps.s.u-tokyo.ac.jp/%7eyokoyama/lecture/suchikaiseki/
「東京大学 横山央明」の検索で横山ホームページ。「講義・演習」リンクからたどれる。
ダウンロードしたファイルを解凍すると、
# ls
advection_ftcs/
advection_lax/
burgers_upwindco/
soundwave/
burgers_upwindnc/
となる。ディレクトリ advection lax については、この解説書の第 5 節以後で詳しく説明する。それ以外のディ
レクトリについては、advection ftcs は実習 1.3 の、burgers upwindnc は実習 2.1 の、burgers upwindnc
は実習 2.2 の、soundwave は実習 3.1 の、それぞれの模範解答プログラムが入っている。参考にしてほしい。
2
I. 実習課題・レポート課題
1
線形波動方程式
実習 1.1
線形波動方程式の差分解法パッケージを、このテキストの第 5 節の記述にそって動かしなさい。サンプル
プログラムは、Lax-Friedrich 法を用いて、１次元線形波動方程式
∂u
∂u
+c
=0
∂t
∂x
(1)
を計算するものである。ただし c = 1 としている。また初期条件は、x ≤ 0.1 に対して u(0, x) = 1、x > 0.1
に対して u(0, x) = 0 である。サンプルプログラムでは、∆x = 0.01、∆t = 0.005、つまり Courant 数は
ν ≡ c∆t/∆x = 0.5。図 1 のような図画が描かれる。なおこの計算に対応する解析解は図 2。
実習 1.2
Lax-Friedrich 法によるプログラムで、時間ステップの跳び ∆t について、CFL 条件を満たす場合と、満た
さない場合とをそれぞれ計算してみなさい。
実習 1.3
このプログラムを FTCS 法によるものに書きかえ、計算結果をグラフ表示しなさい。この結果は図 3 のよ
うになり、数値不安定が起きて解が発散することがわかる。
レポート課題 1.4
このプログラムを 2 段階 Lax-Wendroff 法（人工粘性なし）によるもの、に書きかえ、計算結果をグラフ表
示しなさい。
レポート課題 1.5
このプログラムを 2 段階 Lax-Wendroff 法＋人工粘性によるもの、に書きかえ、計算結果をグラフ表示しな
さい。
2
Burgers 方程式
粘性なしの Burgers 方程式
∂u
∂u
+u
=0
∂t
∂x
を初期条件
u(0, x) = 1
for x < 0.1
u(0, x) = 0 for x ≥ 0.1
のもとで解く。この問題に対しては解析解（弱解）があって、
u(t, x) = 1 for x < t/2 + 0.1
u(t, x) = 0
for x ≥ t/2 + 0.1
3
(2)
2
t=0.0
t=0.1
t=0.2
t=0.3
t=0.4
t=0.5
t=0.6
1.5
1
0.5
0
-0.5
0
0.2
0.4
0.6
0.8
1
図 1: Lax-Friedrich 法による線形波動方程式の計算結果
2
t=0.0
t=0.1
t=0.2
t=0.3
t=0.4
t=0.5
t=0.6
1.5
1
0.5
0
-0.5
0
0.2
0.4
0.6
0.8
1
図 2: 線形波動方程式の解析解
2
t=0.0
t=0.1
t=0.2
t=0.3
t=0.4
t=0.5
t=0.6
1.5
1
0.5
0
-0.5
0
0.2
0.4
0.6
0.8
1
図 3: FTCS 法による線形波動方程式の計算結果。数値不安定の結果、解が発散してプロット範囲の外まで
はみ出している。
となる（図 4）。
実習 2.1
風上差分法を用いて解くプログラムを作成し、結果をグラフ表示しなさい。その際、以下のような、非保存
形の差分で解いてみよ。
∆t n
(u − unj−1 )
∆x j
結果は図 5 のようになり、期待される解析解（図 4）とは大きく異なってしまう。
un+1
= unj − unj
j
4
2
t=0.0
t=0.1
t=0.2
t=0.3
t=0.4
t=0.5
t=0.6
1.5
1
0.5
0
-0.5
0
0.2
0.4
0.6
0.8
1
図 4: 粘性なし Burgers 方程式の解析解
2
t=0.0
t=0.1
t=0.2
t=0.3
t=0.4
t=0.5
t=0.6
1.5
1
0.5
0
-0.5
0
0.2
0.4
0.6
0.8
1
図 5: 粘性なし Burgers 方程式の非保存形を風上差分法で解いた数値解
2
t=0.0
t=0.1
t=0.2
t=0.3
t=0.4
t=0.5
t=0.6
1.5
1
0.5
0
-0.5
0
0.2
0.4
0.6
0.8
1
図 6: 粘性なし Burgers 方程式の保存形を風上差分法で解いた数値解
実習 2.2
同じく風上差分法を用いて解くプログラムを作成し、結果をグラフ表示しなさい。ただし今回は、以下のよ
うな手順で、保存形の差分で解いてみよ。つまり Burgers 方程式の保存形は
( )
∂ u2
∂u
+
=0
∂t
∂x 2
と書けるから、差分式を
un+1
= unj −
j
∆t n
n
(f − fj−1
)
∆x j
fjn =
(unj )2
2
5
(3)
として解いてみよ。結果は図 6 のようになり、期待される解析解（図 4）とだいたい合っている。
レポート課題 2.3
同じ問題を、2 段階 Lax-Wendroff 法＋人工粘性で解きなさい。
3
流体方程式
実習 3.1
等温流体方程式の差分解法パッケージを動かしてみなさい。このサンプルプログラムは、
「2 段階 Lax-Wendroff
法＋人工粘性」で等温流体の方程式
∂ρ
∂
+
(ρVx ) = 0
∂t
∂x
∂
∂
(ρVx ) +
(p + ρVx2 ) = 0
∂t
∂x
(4)
(5)
p = ρRT = ρCs2
(6)
を計算するものである。ただし、温度 T は定数で時間・空間ともに変化しないものとする。つまり p は ρ に
比例し、ここでは比例定数を Cs2 としている。問題設定は、微小振幅擾乱が伝わるようすで、初期条件は、
ρ(x) = 1 + 0.1 cos (2π[x − 0.5]/0.4)
for |x − 0.5| < 0.1
ρ(x) = 1 for |x − 0.5| > 0.1
Vx (x) = 0
for all x
とした。また Cs = 1 としている。結果は図 7 のようになる。
1.2
t=0.0
t=0.1
t=0.2
t=0.3
t=0.4
t=0.5
t=0.6
1.15
1.1
1.05
1
0.95
0.9
0
0.2
0.4
0.6
0.8
1
図 7: 等温流体方程式の数値解
レポート課題 3.2
等温流体の方程式のプログラムを参考にして、（断熱）流体の方程式
∂
∂ρ
+
(ρVx ) = 0
∂t
∂x
(7)
∂
∂
(ρVx ) +
(p + ρVx2 ) = 0
∂t
∂x
(8)
6
∂
∂ϵ
+
[(ϵ + p)Vx ] = 0
∂t ∂x
{
}
1 2
p = (γ − 1) ϵ − ρVx
2
(9)
(10)
を「2 段階 Lax-Wendroff 法＋人工粘性」を用いて解くプログラムを作成しなさい。これを用いて以下の、
衝撃波管問題（「Sod の問題」ともいう）を解きなさい。ただし比熱比 γ = 1.4 とする。
（この問題に対して
も解析解が存在する。「レポート挑戦課題 3.3」参照。）
0 < x < 0.5 に対して、ρ = 1、p = 1、Vx = 0。
0.5 < x < 1 に対して、ρ = 0.125、p = 0.1、Vx = 0。
t = 0.14154 のときの ρ、p、Vx の分布を出力すること。
レポート挑戦課題 3.3
［この課題は、成績評価には含めない。興味あるひとは、挑戦してみてください。］課題 3.2 の問題に対す
る解析解を求めよ。解法については、松元亮治ほか「天体とスペースプラズマのシミュレーションサマース
クール講義テキスト」http://www.astro.phys.s.chiba-u.ac.jp/netlab/summer-school 2004/textbook.html
の第 1.2 節「システム方程式の解法」（紙出力版では第 2 章）で詳しく解説されている。
4
補足解説 — 保存形による差分法
保存形で書かれた方程式
∂u ∂f
+
=0
∂t
∂x
を差分化すると
un+1
= unj −
j
∆t ˜n+1/2 ˜n+1/2
(f
− fj−1/2 )
∆x j+1/2
(11)
(12)
n+1/2
となる。f˜j+1/2 は数値流束である。
FTCS 法: 数値流束は
1 n
n+1/2
f˜j+1/2 = (fj+1
+ fjn )
2
(13)
fjn = f (unj )
(14)
[
]
∆x n
1
n+1/2
n
f˜j+1/2 =
(fj+1
+ fjn ) −
(uj+1 − unj )
2
∆t
(15)
ただし、
である（以下この節では同様）。
Lax-Friedrich 法: 数値流束は
となる。
2 段階 Lax-Wendroff 法: 数値流束は
n+1/2
uj+1/2 =
1 ∆t n
1 n
(uj+1 + unj ) −
(f
− fjn )
2
2 ∆x j+1
n+1/2
n+1/2
f˜j+1/2 = f (uj+1/2 )
un+1
= unj −
j
∆t ˜n+1/2 ˜n+1/2
(f
− fj−1/2 )
∆x j+1/2
となる。
7
(16)
(17)
(18)
II. 配布パッケージの説明
配布パッケージの動かし方
5
線形波動方程式（移流方程式 advection equation）の差分解法を試してみましょう。作業はディレクトリ
advection lax で行ないます。ここには以下のファイルが用意されています。
# cd
# tar xvf jisshu-sample.tar
# cd jisshu-sample
# ls
advection_ftcs
burgers_upwindco
soundwave
advection_lax
burgers_upwindnc
# cd advection_lax
# ls
main.f
plot.gpl
プログラムは Fortran 言語を用いて書かれています。上記のリストで main.f が Fortran プログラムファイ
ルです。plot.gpl は可視化に使うファイルです。
5.1
プログラムのコンパイルと実行
プログラムをコンパイル（機械語翻訳）する方法について説明します。コマンド gfortran を実行します。
するとプログラムがコンパイルされます。
# gfortran main.f
# ls
a.out
main.f
plot.gpl
正しくコンパイルされると実行ファイル a.out を作成します。
以下のようにして機械語翻訳されたプログラム（ファイル a.out）を実行します。
# ./a.out
write
write
....
step=
step=
0 time= 0.000E+00 nd =
20 time= 0.100E+00 nd =
1
2
write
stop
step=
step=
121 time= 0.605E+00 nd =
121 time= 0.605E+00
8
### normal stop ###
# ls
a.out
main.f
out.dat
plot.gpl
正しく実行されるとデータファイル out.dat を出力します。
8
出力ファイルの説明 (out.dat)
5.2
出力ファイル out.dat はアスキー形式（人間が読める文字形式）でかかれており、ファイルを直接エディ
タで開いて見ることができます。ファイルをみてみましょう。第 1 行目から第 100 行目に渡って、始めの
時間データの x 座標 (x)、そこでの変数の値 (u) が順にかかれています。第 101 行目以降に、次の時間デー
タが書かれています。書式については、Fortran プログラム main.f に Format 文で指定されていますので、
参考にしてください。
# head out.dat
0.000E+00 0.100E+01
0.100E-01 0.100E+01
0.200E-01
0.300E-01
0.100E+01
0.100E+01
0.400E-01
0.500E-01
(後略)
0.100E+01
0.100E+01
5.3
結果の可視化表示
結果の表示には gnuplot 可視化プログラムを利用します。
5.3.1
gnuplot の起動 (gnuplot)
まずは gnuplot を実行してみましょう。
# gnuplot
すると以下のようになり、gnuplot が起動します。
G N U P L O T
Version ...
....
gnuplot>
5.3.2
データ表示 (plot.gpl)
データの表示には plot.gpl というファイルに書かれたプログラムを用います。以下のように入力してみて
下さい。
gnuplot> load ’plot.gpl’
図 1 のような図画が描かれます。なおこの計算に対応する解析解は図 2 です。
9
5.3.3
gnuplot の終了 (quit)
quit を入力すると gnuplot を終了することができます。
gnuplot> quit
プログラムの変更について
6
6.1
時間ステップの跳び (dt)
時間ステップの跳びは、89 行目の dt の値を変更します。この値は、数値スキームの安定性条件（CFL 条
件）を満たすように決めなければなりません。
6.2
dt=0.5e-2
メッシュ数の設定 (ix) 空間座標の跳び (dx)
メッシュ数を変更するには、6 行目の parameter 文にある ix の値を変更します。メッシュ数をかえること
で数値計算の分解能をあげることができます。ただし、CFL 条件をやぶらないように、同時に dt も調整す
る必要があることに注意してください。
parameter (ix=100)
また空間座標の跳びを変更するには 54 行目の以下の箇所を変更します。このプログラムでは、計算領域全
体の長さが「1」になるように設定してあります。
6.3
dx=1.e0/real(ix)
最終ステップ数、出力の設定 (tend, dtout)
最終時刻と出力データの間隔は、32 行目の tend と 33 行目の dtout の値をそれぞれ変更します。
c
time control parameters
tend=0.6e0
6.4
dtout=0.1e0
計算エンジンの変更
サンプルプログラムの計算エンジンは Lax-Friedrich 法になっています。96 行目あたりにある下記の部分
で、計算エンジンを変更してください。
10
c Lax-Friedrich
do i=1,ix
f0(i)=cs*u(i)
enddo
do i=1,ix-1
f(i)=0.5e0*(f0(i+1)+f0(i))-0.5e0*dx/dt*(u(i+1)-u(i))
enddo
do i=2,ix-1
um(i)=u(i)-dt/dx*(f(i)-f(i-1))
enddo
do i=2,ix-1
u(i)=um(i)
enddo
11
プログラムの中身
7
サンプルプログラム、main.f
c======================================================================|
c
array definitions
c======================================================================|
implicit none
integer ix,i
parameter (ix=100)
real x,dx
dimension x(ix)
real u,f,um,f0
dimension u(ix),f(ix),um(ix),f0(ix)
real time,dt,tend
integer ns,nstop
real tout,dtout
integer nd
real cs
c======================================================================|
c
prologue
c======================================================================|
c----------------------------------------------------------------------|
c
initialize
do i=1,ix
x(i)=0.e0
u(i)=0.e0
f(i)=0.e0
um(i)=0.e0
f0(i)=0.e0
enddo
c----------------------------------------------------------------------|
c
time control parameters
tend=0.6e0
dtout=0.1e0
nstop=1000000
12
c----------------------------------------------------------------------|
c file open
open(10,file=’out.dat’,form=’formatted’)
c----------------------------------------------------------------------|
c initialize counters
nd=1
time = 0.e0
tout = 0.e0
ns=0
c----------------------------------------------------------------------|
c
setup numerical model (grid, initial conditions, etc.)
cs=1.e0
dx=1.e0/real(ix)
x(1)=0.e0
do i=2,ix
x(i)=x(i-1)+dx
enddo
do i=1,ix
if (x(i).le.0.1e0) then
u(i)=1.e0
else
u(i)=0.e0
endif
enddo
c----------------------------------------------------------------------|
c
data output
do i=1,ix
write(10,910) x(i),u(i)
enddo
tout=tout+dtout
write(6,913) ns,time,nd
nd=nd+1
13
c======================================================================|
c
time integration
c======================================================================|
1000 continue
ns = ns+1
c----------------------------------------------------------------------|
c
time spacing
dt=0.005e0
time = time+dt
c----------------------------------------------------------------------|
c
solve difference equations
c Lax-Friedrich
do i=1,ix
f0(i)=cs*u(i)
enddo
do i=1,ix-1
f(i)=0.5e0*(f0(i+1)+f0(i))-0.5e0*dx/dt*(u(i+1)-u(i))
enddo
do i=2,ix-1
um(i)=u(i)-dt/dx*(f(i)-f(i-1))
enddo
do i=2,ix-1
u(i)=um(i)
enddo
c----------------------------------------------------------------------|
c
boundary condition
u(1)=u(2)
u(ix)=u(ix-1)
14
c----------------------------------------------------------------------|
c
data output
if (time.ge.tout) then
do i=1,ix
write(10,910) x(i),u(i)
enddo
tout=tout+dtout
write(6,913) ns,time,nd
nd=nd+1
endif
c----------------------------------------------------------------------|
c
loop test
if (ns .lt. nstop .and. time .lt. tend) goto 1000
c======================================================================|
c
epilogue
c======================================================================|
c----------------------------------------------------------------------|
c
data output
do i=1,ix
write(10,910) x(i),u(i)
enddo
write(6,913) ns,time,nd
c----------------------------------------------------------------------|
c
file close
close(10)
c----------------------------------------------------------------------|
c
ending message
write(6,915) ns,time
write(6,*) ’ ### normal stop ###’
stop
15
910
913
915
format (1x,e10.3,1x,e10.3)
format (1x,’ write
’,’step=’,i8,’ time=’,e10.3,’ nd =’,i3)
format (1x,’ stop
’,’step=’,i8,’ time=’,e10.3)
end
サンプルプログラム、plot.gpl
# if plot to PostScript file then set to 1 othewise set to 0
plot_ps=0
# for Postscript plot
if (plot_ps == 1) \
set term postscript portrait 18; \
set output "gnuplot.ps"; \
set size 1.0,0.5
set xrange [0:1]
set yrange [-0.5:2]
plot 0 linetype 3 notitle, \
"out.dat" every ::0::99 using 1:2 with lines linewidth 3 title "t=0.0", \
"out.dat" every ::100::199 using 1:2 with lines linewidth 3 title "t=0.1", \
"out.dat" every ::200::299 using 1:2 with lines linewidth 3 title "t=0.2", \
"out.dat" every ::300::399 using 1:2 with lines linewidth 3 title "t=0.3", \
"out.dat" every ::400::499 using 1:2 with lines linewidth 3 title "t=0.4", \
"out.dat" every ::500::599 using 1:2 with lines linewidth 3 title "t=0.5", \
"out.dat" every ::600::699 using 1:2 with lines linewidth 3 title "t=0.6"
# for Postscript plot
# set string "x11" to "aqua" for MacOS
if (plot_ps == 1) \
set terminal x11; \
set size 1.0,1.0
16

コンピュータ物理学2 第1回 (2015.10.2)

スライド 1

グラフのPLOT方法 - Oyanagi Laboratory WWW Server

TEL：06-6718-5211 FAX：06-6251-3308 CD マスタリング指示書

gnuplot でアニメーション - IT pass (Informational

並列計算機論 - 奈良女子大学理学部情報科学科

情報処理概論Ⅰ 2007 第3回

PDF ファイルへの 3D モデルの埋め込み

スライド 1

LU分解法（2）

高次保存量を持つCAの確率化について

Gnuplotとの連携

課題 2: 自己組織化モデル

予稿集

Application Note 資料作成 040520 資料番号 68K - accelmpc.co.jp

LU分解法（1）

石油タンクの地震被害想定用エキスパートシステム の作成に関する研究

06/01 統計用言語 R の使い方

ルンゲクッタ型多項式の安定性

引張軟化曲線多直線近似解析 実行プログラム（exe ﾌｧｲﾙ）概要

Preconditioned Iterative Linear Solvers for

並列化による高速化技法