葛西 真寿 弘前大学大学院理工学研究科
gnuplot は，基本はグラフを作成するためのアプリケーションソフトウェアですが，最近のバージョンでは，繰り返し処理や三項演算子を使った再帰定義を利用することで，簡単なプログラミングによる数値解析ができます。プログラミングに必要なコマンドと，いくつかの応用例ををまとめておきます。
書式
if (条件式) {
条件式が満たされる場合に実行する文
複数行も可
} else {
条件式が満たされない場合に実行する文
複数行も可
}x = rand(0) # 開区間 (0:1) 内の疑似乱数値を返す
if (x <= 0.5) {
print x, " は 0.5 以下"
} else {
print x, " は 0.5 より大きい"
}
$\displaystyle \sum_{i=1}^{5} i$ を計算する例
sum = 0
do for [i = 1: 5] {
sum = sum + i
print "1 から ", i, " までの和は ", sum
}
sum = 0
i = 1
while (i<=5) {
sum = sum + i
print "1 から ", i, " までの和は ", sum
i = i + 1
}
$\mbox{isum}(n) = \displaystyle \sum_{i=1}^{n} i$ を以下のように考えます。
もし，$n = 1$ なら $\mbox{isum}(n) = \mbox{isum}(1) = 1$
もし，$n > 1$ なら $\mbox{isum}(n) = \mbox{isum}(n-1) + n$
これを gnuplot の「三項演算子」を使って以下のように定義します。
# 三項演算子を使った再帰定義
# isum(n) が自身の1つ前の値 isum(n-1) を使って定義されています
isum(n) = (n == 1 ? 1 : isum(n-1) + n)
# i = 1 から 5 まで増分 2 ごとに print 文を実行する例
do for [i = 1: 5: 2] {
print "1 から ", i, " までの和は ", isum(i)
}
解析的な微分の定義は（前方差分の極限として） $$\frac{df}{dx} \equiv \lim_{h \rightarrow 0} \frac{f(x+h) - f(x)}{h}$$ でした。滑らかな関数であれば，（後方差分の極限として） $$\frac{df}{dx} = \lim_{h \rightarrow 0} \frac{f(x) - f(x-h)}{h}$$ としても良いです。
現実世界の gnuplot では，$ h\rightarrow 0$ の極限はとれませんから十分小さい値として h を定義し，近似的な数値微分を（前方差分と後方差分の平均である中心差分として）以下のように定義しましょう。
以下の例では，$f(x) = \sin x$ に対して $f'(x) = \cos x$ と，上記の数値微分との誤差を描いています。
f(x) = sin(x)
# f(x) 微分
diff_f(x, h) = (f(x+h)-f(x-h))/(2.0*h)
set samples 1000
plot [-2*pi:2*pi] cos(x) - diff_f(x, 1.0E-6)
plot [-2*pi:2*pi] cos(x) - diff_f(x, 1.0E-7)
シンプソン法で $\displaystyle \int_a^b f(x)\, dx$ を求める。
積分区間 $[a, b]$ を $N (=2n)$ 等分（偶数等分）し， $$h = \frac{b-a}{N}, \quad f_i = f(a + i\,h)$$ とする。シンプソン法の公式は \begin{eqnarray} \int_a^b f(x)\, dx &\simeq& \frac{h}{3}( f_0 + 4 f_1 + f_2) + \frac{h}{3}( f_2 + 4 f_3 + f_4) + \cdots\\ &&\quad \cdots + \frac{h}{3}( f_{N-2} + 4 f_{N-1} + f_N)\\ &=&\frac{h}{3} \left(f_0 + f_{2n} + 2 \sum_{i=1}^{n-1} f_{2i} + 4 \sum_{i=1}^{n} f_{2n-1}\right) \end{eqnarray}
$\displaystyle \int_0^{\pi/2} \sin x \,dx$ の定積分を，積分区間を $N=20$ 分割してシンプソン法で解く例です。
do for 文の繰り返し処理で和をとっています。
f(x) = sin(x) # 被積分関数
a = 0.0 # 積分の下端。実数で入力。
b = pi/2 # 積分の上端。実数で入力。
N = 10 # 分割数 N は偶数とすること（シンプソン法の決め事）
# ====
h = (b-a)*1./N # * 1. は実数で入力しない人対策
xi(i) = a + i*h
fi(i) = f(xi(i))
simpson_f = 0
do for [i=2: N: 2] {
simpson_f = simpson_f + h/3.*(fi(i-2) + 4.*fi(i-1) + fi(i))
}
print N," 分割の数値解は ", simpson_f
# ==== 精度確認のため，分割数を倍にして計算
N = 2*N
h = (b-a)*1./N
xi(i) = a + i*h
fi(i) = f(xi(i))
simpson_f = 0
do for [i=2: N: 2] {
simpson_f = simpson_f + h/3.*(fi(i-2) + 4.*fi(i-1) + fi(i))
}
print N," 分割の数値解は ", simpson_f
$\displaystyle \int_0^{\pi/2} \sin x \,dx$ の定積分を，積分区間を $N=20$ 分割してシンプソン法で解く例です。
再帰的定義関数で和をとっています。
以下の simp_f() の定義文のように，文が長くなって1行で収まらない場合は，継続を表す \ で改行できます。
undefine a b N i
f(x) = sin(x) # 被積分関数
# a 積分の下端。実数で入力。
# b 積分の上端。実数で入力。
# N 分割数 N は偶数とすること（シンプソン法の決め事）
simp_f(a, b, N, i) = (i == 2 ? \
(b-a)/(N*3.0)*(f(a+(i-2)*(b-a)/N) + 4.0*f(a+(i-1)*(b-a)/N) + f(a+i*(b-a)/N)) : \
simp_f(a, b, N, i-2) + \
(b-a)/(N*3.0)*(f(a+(i-2)*(b-a)/N) + 4.0*f(a+(i-1)*(b-a)/N) + f(a+i*(b-a)/N)))
# 実際の計算の際は，simp_f(a, b, N, N) と同じ N を2回書く
N = 10
print N," 分割の数値解は ", simp_f(0.0, pi/2, N, N)
N = 20
print N," 分割の数値解は ", simp_f(0.0, pi/2, N, N)
$f(x) = 0$ の解を求める。
まず，何らかの方法で（例えば，plot f(x) としてグラフを描いてみて $x$ 軸と交差する点のあたりをつけるなどして）$f(x) = 0$ の解である $x$ に近いと思われる値 $x_k$ がわかったとする。
この $x_k$ のまわりに $f(x)$ をテイラー展開すると， $$ f(x) \simeq f(x_k) + f'(x_k)\,(x - x_k)$$
$f(x) = 0$ であるから，上式の左辺をゼロとおいて，$x$ について解くと， $$ x \equiv x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}$$
これを $|f(x_{k+1})| < \epsilon$ （ここで $\epsilon$ 十分小さい数値，たとえば $1.0\times 10^{-8}$ とか？）になるまで反復計算を行い，このときの $x_{k+1}$ をもって，$f(x) = 0$ の近似解とする。
$f'(x_k)$ を求めるのが面倒な場合は，以下のようにして数値的な微分で代用する。 $$f'(x_k) \simeq \frac{f(x_k + h) - f(x_k - h)}{2 h}, \quad |h| \ll 1$$
つまり， $$x_{k+1} = x_k - \frac{f(x_k)\times 2 h}{f(x_k + h) - f(x_k - h)}$$
まず，解を求める関数 $f(x)$ を定義し，$f(x) = 0$ となる $x$ のあたりをつけるためにグラフを描いてみます。
以下では例として組み込み関数の1つである besj0(x) のゼロ点を求めています。gnuplot の組み込み関数については，たとえば以下を参照。
f(x) = besj0(x) # 解（ゼロ点）を求める関数
set grid
set zeroaxis
plot [0:10] f(x) title "J_0 ベッセル関数"
$f(x) = 0$ となる $x$ の1つは，$2$ から $3$ の間にあることがわかります。
探索の初期値である x は上図から推定される 2.5 を入れてみます。解の収束条件 h は，数値解の精度に関係しますので，いくつか値を変えて出力させます。
f(x) = besj0(x)
# x には探索の初期値を入れる
# h は解の収束条件である微小量と数値微分をするときに使う微小量を兼ねる。
x = 2.5
h = 1E-6
while (abs(f(x)) > h) {
x = x - f(x)*2.*h/(f(x + h)-f(x - h))
}
xx1 = x
print xx1
# ==== 精度確認のため，h を変えて計算
x = 2.5
h = 1E-12
while (abs(f(x)) > h) {
x = x - f(x)*2.*h/(f(x + h)-f(x - h))
}
xx2 = x
print xx1
print xx2
以下のように，再帰定義により $f(x) = 0$ の解を求める関数 newton_f(x, h) を定義します。
探索の初期値である x は上図から推定される 2.5 を入れてみます。解の収束条件 h は，数値解の精度に関係しますので，いくつか値を変えて出力させます。
# 再帰定義による解
# x には探索の初期値を入れる
# h は解の収束条件である微小量と数値微分をするときに使う微小量を兼ねる。
newton_f(x, h) = \
(abs(f(x)) > h ? newton_f(x-f(x)*2.*h/(f(x+h)-f(x-h)), h) : x)
# 解を表示する
print newton_f(2.5, 1E-6)
print newton_f(2.5, 1E-8)
print newton_f(2.5, 1E-10)
print newton_f(2.5, 1E-12)
print newton_f(2.5, 1E-14)
解の収束条件 h を変えても変わらない桁までが，精度のよい数値解ということになります。上の例では，
2.404825557695
くらいでしょうか。
print besj0(2.404825557695)
4次のルンゲ・クッタ法で常微分方程式を解きます。
1階の常微分方程式 $$ \frac{dx}{dt} = f(x, t)$$ を4次のルンゲ・クッタ法で解く例です。
初期条件を $t=t_0$ で $x(t_0) = x_0$ とし，$t = t_1$ までを $N$ 分割して解きます。
$$ h = \frac{t_1 - t_0}{N}, \quad t_i = t_0 + i\, h, \quad x_i = x(t_i)$$とし，以下の式から次のステップの値 $x_{i+1}$ を決めます。
\begin{eqnarray} k_1 &=& h\, f(x_i, t_i) \\ k_2 &=& h\, f(x_i + k_1/2, t_i + h/2) \\ k_3 &=& h\, f(x_i + k_2/2, t_i + h/2) \\ k_4 &=& h\, f(x_i + k_3, t_i + h) \\ x_{i+1} &=& x_i + \frac{1}{6} (k_1 + 2 k_2 + 2 k_3 + k_4) \end{eqnarray}簡単なロジスティック方程式 $$ \frac{dx}{dt} = (1-x)\,x$$ を，初期条件 $t_0 = 0$ で $x_0 = x(t_0) = 0.1$，$t = t_0 = 0$ から $t_1 = 10$ までを $ N = 100$ 分割して解きます。
f(x, t) = (1-x)*x
t0 = 0.0
x0 = 0.1
t1 = 10.0
N = 100
h = (t1 - t0) * 1.0/N
# 計算結果を配列に入れます。
array T[N]
array X[N]
t = t0; x = x0
do for [i=1: N]{
k1 = h*f(x, t)
k2 = h*f(x+0.5*k1, t+0.5*h)
k3 = h*f(x+0.5*k2, t+0.5*h)
k4 = h*f(x+k3, t+h)
X[i] = x + (k1 + 2.*k2 + 2.*k3 + k4)/6
T[i] = t + h
x = X[i]
t = T[i]
}
配列に入れた計算結果をファイル rk.dat に書き出し，そのファイルを読み込んで plot する例です。
set print "rk.dat" # ファイル名の指定。以後はこのファイルに print
print "# t x"
print t0, " ", x0
do for [i=1:N]{
print T[i], " ", X[i]
}
set print # 書き終わったら標準出力に戻す。
# ファイル "rk.dat" の数値データをプロットします。
set yrange [0:1]
set pointsize 0.3 # 点の大きさを小さめに
plot "rk.dat" pt 7 lc 7 # pointtype と 色 (lc) を適宜設定
plot X using (T[$1]):2 pt 7 notitle
分割数 N を変えて結果をみます。分割数を変えても変わらない桁数までをもって，精度の良い数値解とします。
f(x, t) = (1-x)*x
t0 = 0.0
x0 = 0.1
t1 = 10.0
do for [i=2:8:2]{
N = 100*i
h = (t1 - t0) * 1.0/N
# 計算結果を配列に入れます。
array T[N]
array X[N]
t = t0; x = x0
do for [i=1: N]{
k1 = h*f(x, t)
k2 = h*f(x+0.5*k1, t+0.5*h)
k3 = h*f(x+0.5*k2, t+0.5*h)
k4 = h*f(x+k3, t+h)
X[i] = x + (k1 + 2.*k2 + 2.*k3 + k4)/6
T[i] = t + h
x = X[i]
t = T[i]
}
print N, " 分割時の数値 ", T[N], " ", X[N]
}
数値解析の本来の目的は，解析的に解けない問題を何とかして解きたい，解いたときの精度もちゃんと把握したい，ということです。
なので，上記の数値解がどの程度の精度で精確であるかを調べるために解析解と比較するのは本来は邪道なのですが，たまたま，上記の微分方程式は解析的に解けて，以下のようになります。
$$x(t) = \frac{\exp(t)}{9 + \exp(t)}$$次のような2階常微分方程式をルンゲ・クッタ法で解く。 $$ \frac{d^2 x}{dt^2 } = f\left(x, \frac{dx}{dt}, t\right)$$
この場合には， $\displaystyle v \equiv \frac{dx}{dt}$ とおけば，次のような連立1階常微分方程式の形に帰着できる。 \begin{eqnarray} \frac{dx}{dt} &=& F_1(x, v, t) = v \\ \frac{dv}{dt} &=& F_2(x, v, t) = f(x, v, t) \end{eqnarray}
初期条件を $t=t_0$ で $x(t_0) = x_0, v(t_0) = v_0$ とし，$t = t_1$ までを $N$ 分割して解きます。
$$ h = \frac{t_1 - t_0}{N}, \quad t_i = t_0 + i\, h, \quad x_i = x(t_i), \quad v_i = v(t_i)$$とし，以下の式から次のステップの値 $x_{i+1}, v_{i+1}$ を決めます。
\begin{eqnarray} k_1 &=& h\, F_1(x_i, v_i, t_i) \\ m_1 &=& h\, F_2(x_i, v_i, t_i) \\ k_2 &=& h\, F_1(x_i + k_1/2, v_i + m_1/2, t_i + h/2) \\ m_2 &=& h\, F_2(x_i + k_1/2, v_i + m_1/2, t_i + h/2) \\ k_3 &=& h\, F_1(x_i + k_2/2, v_i + m_2/2, t_i + h/2) \\ m_3 &=& h\, F_2(x_i + k_2/2, v_i + m_2/2, t_i + h/2) \\ k_4 &=& h\, F_1(x_i + k_3, v_i + m_3, t_i + h) \\ m_4 &=& h\, F_2(x_i + k_3, v_i + m_3, t_i + h) \\ x_{i+1} &=& x_i + \frac{1}{6} (k_1 + 2 k_2 + 2 k_3 + k_4)\\ v_{i+1} &=& v_i + \frac{1}{6} (m_1 + 2 m_2 + 2 m_3 + m_4) \end{eqnarray}簡単な減衰＋強制振動の方程式
$$ \frac{d^2 x}{dt^2} = -x - a \frac{dx}{dt} + b \cos(t) $$を，
して解きます。
$N$ の値をいくつか変えてみて，精度を確認します。
a = 0.5
b = 0.2
F1(x, v, t) = v
F2(x, v, t) = -x - a*v + b*cos(t)
t0 = 0
x0 = 3.0
v0 = 0.0
t1 = 20
do for [j=0: 3]{
N = 100*2**j
h = (t1 - t0) * 1.0/N
# 計算結果を配列に入れます。
array T[N]
array X[N]
array V[N]
t = t0; x = x0; v = v0
do for [i=1: N]{
k1 = h*F1(x, v, t)
m1 = h*F2(x, v, t)
k2 = h*F1(x+0.5*k1, v+0.5*m1, t+0.5*h)
m2 = h*F2(x+0.5*k1, v+0.5*m1, t+0.5*h)
k3 = h*F1(x+0.5*k2, v+0.5*m2, t+0.5*h)
m3 = h*F2(x+0.5*k2, v+0.5*m2, t+0.5*h)
k4 = h*F1(x+k3, v+m3, t+h)
m4 = h*F2(x+k3, v+m3, t+h)
X[i] = x + (k1 + 2.*k2 + 2.*k3 + k4)/6
V[i] = v + (m1 + 2.*m2 + 2.*m3 + m4)/6
T[i] = t + h
t = T[i]
x = X[i]
v = V[i]
}
print X[N]," N = ", N
}
$N = 400$ で小数点以下7桁くらいの精度が出ているようですので，$N=400$ として続けます。
まずは，$a = 0, b = 0$，すなわち単振動の場合の方程式を解き，結果をファイルに保存します。
a = 0
b = 0
F1(x, v, t) = v
F2(x, v, t) = -x - a*v + b*cos(t)
t0 = 0
x0 = 3.0
v0 = 0.0
t1 = 20
N = 400
h = (t1 - t0) * 1.0/N
# 計算結果を配列に入れます。
array T[N]
array X[N]
array V[N]
t = t0; x = x0; v = v0
do for [i=1: N]{
k1 = h*F1(x, v, t)
m1 = h*F2(x, v, t)
k2 = h*F1(x+0.5*k1, v+0.5*m1, t+0.5*h)
m2 = h*F2(x+0.5*k1, v+0.5*m1, t+0.5*h)
k3 = h*F1(x+0.5*k2, v+0.5*m2, t+0.5*h)
m3 = h*F2(x+0.5*k2, v+0.5*m2, t+0.5*h)
k4 = h*F1(x+k3, v+m3, t+h)
m4 = h*F2(x+k3, v+m3, t+h)
X[i] = x + (k1 + 2.*k2 + 2.*k3 + k4)/6
V[i] = v + (m1 + 2.*m2 + 2.*m3 + m4)/6
T[i] = t + h
t = T[i]
x = X[i]
v = V[i]
}
計算結果をファイルに書き出します。
set print "rk2-a00-b00.dat"
print "# a = ", a, " b = ", b
print "# t x v"
do for [i=1:N]{
print T[i], " ", X[i], " ", V[i]
}
set print
! head rk2-a00-b00.dat
ファイル rk2-a00-b00.dat のデータを plot します。
set grid
set zeroaxis
plot "rk2-a00-b00.dat" w l title "a=0, b=0"
次の例として，$a = 0.5, b = 0.2$ の場合を解き，結果をファイルに数値データとして保存します。
a = 0.5
b = 0.2
F1(x, v, t) = v
F2(x, v, t) = -x - a*v + b*cos(t)
t0 = 0
x0 = 3.0
v0 = 0.0
t1 = 20
N = 400
h = (t1 - t0) * 1.0/N
# 計算結果を配列に入れます。
array T[N]
array X[N]
array V[N]
t = t0; x = x0; v = v0
do for [i=1: N]{
k1 = h*F1(x, v, t)
m1 = h*F2(x, v, t)
k2 = h*F1(x+0.5*k1, v+0.5*m1, t+0.5*h)
m2 = h*F2(x+0.5*k1, v+0.5*m1, t+0.5*h)
k3 = h*F1(x+0.5*k2, v+0.5*m2, t+0.5*h)
m3 = h*F2(x+0.5*k2, v+0.5*m2, t+0.5*h)
k4 = h*F1(x+k3, v+m3, t+h)
m4 = h*F2(x+k3, v+m3, t+h)
X[i] = x + (k1 + 2.*k2 + 2.*k3 + k4)/6
V[i] = v + (m1 + 2.*m2 + 2.*m3 + m4)/6
T[i] = t + h
t = T[i]
x = X[i]
v = V[i]
}
set print "rk2-a05-b02.dat"
print "# a = ", a, " b = ", b
print "# t x v"
do for [i=1:N]{
print T[i], " ", X[i], " ", V[i]
}
set print
2つの数機計算の結果を重ねてグラフ化します。
replot "rk2-a05-b02.dat" w l title "a=0.5, b=0.2"