はじめに:gnuplot で数値解析

葛西 真寿 弘前大学大学院理工学研究科

gnuplot は,基本はグラフを作成するためのアプリケーションソフトウェアですが,最近のバージョンでは,繰り返し処理や三項演算子を使った再帰定義を利用することで,簡単なプログラミングによる数値解析ができます。プログラミングに必要なコマンドと,いくつかの応用例ををまとめておきます。

条件分岐

書式

if (条件式) {
    条件式が満たされる場合に実行する文
    複数行も可
  } else {
    条件式が満たされない場合に実行する文
    複数行も可
}
In [1]:
x = rand(0)  # 開区間 (0:1) 内の疑似乱数値を返す
if (x <= 0.5) {
    print x, " は 0.5 以下"
  } else {
    print x, " は 0.5 より大きい"
}
Out[1]:
0.222457440974512 は 0.5 以下

繰り返し処理

$\displaystyle \sum_{i=1}^{5} i$ を計算する例

do for

In [2]:
sum = 0
do for [i = 1: 5] {
  sum = sum + i
  print "1 から ", i, " までの和は   ", sum
}
Out[2]:
1 から 1 までの和は   1
1 から 2 までの和は   3
1 から 3 までの和は   6
1 から 4 までの和は   10
1 から 5 までの和は   15

while

In [3]:
sum = 0
i = 1
while (i<=5) {
  sum = sum + i
  print "1 から ", i, " までの和は   ", sum
  i = i + 1
}
Out[3]:
1 から 1 までの和は   1
1 から 2 までの和は   3
1 から 3 までの和は   6
1 から 4 までの和は   10
1 から 5 までの和は   15

三項演算子と再帰定義関数

$\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 の「三項演算子」を使って以下のように定義します。

In [4]:
# 三項演算子を使った再帰定義
# 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)
}
Out[4]:
1 から 1 までの和は   1
1 から 3 までの和は   6
1 から 5 までの和は   15

数値微分

解析的な微分の定義は(前方差分の極限として) $$\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 を定義し,近似的な数値微分を(前方差分と後方差分の平均である中心差分として)以下のように定義しましょう。

\begin{eqnarray} \frac{df}{dx} &\simeq& \frac{1}{2} \left\{\frac{f(x+h) - f(x)}{h} + \frac{f(x) - f(x-h)}{h} \right\} \\ &=& \frac{f(x+h) - f(x-h)}{2h} \end{eqnarray}

以下の例では,$f(x) = \sin x$ に対して $f'(x) = \cos x$ と,上記の数値微分との誤差を描いています。

In [5]:
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)
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 -2x10-10 -1.5x10-10 -1x10-10 -5x10-11 0 5x10-11 1x10-10 1.5x10-10 2x10-10 -6 -4 -2 0 2 4 6 cos(x) - diff_f(x, 1.0E-6) cos(x) - diff_f(x, 1.0E-6)
In [6]:
plot [-2*pi:2*pi] cos(x) - diff_f(x, 1.0E-7)
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 -3x10-9 -2x10-9 -1x10-9 0 1x10-9 2x10-9 3x10-9 -6 -4 -2 0 2 4 6 cos(x) - diff_f(x, 1.0E-7) 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}

do for 文による解法

$\displaystyle \int_0^{\pi/2} \sin x \,dx$ の定積分を,積分区間を $N=20$ 分割してシンプソン法で解く例です。

do for 文の繰り返し処理で和をとっています。

In [7]:
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
Out[7]:
10 分割の数値解は 1.0000033922209
20 分割の数値解は 1.00000021154659

再帰定義関数による解法

$\displaystyle \int_0^{\pi/2} \sin x \,dx$ の定積分を,積分区間を $N=20$ 分割してシンプソン法で解く例です。

再帰的定義関数で和をとっています。

以下の simp_f() の定義文のように,文が長くなって1行で収まらない場合は,継続を表す \ で改行できます。

In [8]:
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)
Out[8]:
10 分割の数値解は 1.0000033922209
20 分割の数値解は 1.00000021154659

方程式の数値解

$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 の組み込み関数については,たとえば以下を参照。

In [9]:
f(x) = besj0(x) # 解(ゼロ点)を求める関数

set grid
set zeroaxis
plot [0:10] f(x) title "J_0 ベッセル関数"
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 0 2 4 6 8 10 J0 ベッセル関数 J0 ベッセル関数

while 文による解法

$f(x) = 0$ となる $x$ の1つは,$2$ から $3$ の間にあることがわかります。

探索の初期値である x は上図から推定される 2.5 を入れてみます。解の収束条件 h は,数値解の精度に関係しますので,いくつか値を変えて出力させます。

In [10]:
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
Out[10]:
2.40482459179177
In [11]:
# ==== 精度確認のため,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
Out[11]:
2.40482459179177
2.40482555769576

再帰定義関数による解法

以下のように,再帰定義により $f(x) = 0$ の解を求める関数 newton_f(x, h) を定義します。

探索の初期値である x は上図から推定される 2.5 を入れてみます。解の収束条件 h は,数値解の精度に関係しますので,いくつか値を変えて出力させます。

In [12]:
# 再帰定義による解
# 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)
Out[12]:
2.40482459179177
2.40482555769559
2.40482555769544
2.40482555769576
2.40482555769578

解の収束条件 h を変えても変わらない桁までが,精度のよい数値解ということになります。上の例では,

2.404825557695

くらいでしょうか。

In [13]:
print besj0(2.404825557695)
Out[13]:
4.01086468266012e-13

微分方程式の数値解

4次のルンゲ・クッタ法で常微分方程式を解きます。

1階常微分方程式

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$ 分割して解きます。

In [14]:
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]
}

配列データの書き出しと plot

配列に入れた計算結果をファイル rk.dat に書き出し,そのファイルを読み込んで plot する例です。

In [15]:
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) を適宜設定
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 7 8 9 10 "rk.dat" "rk.dat"

配列データを直接 plot

配列データを直接 plot する例です。

T[i] を横軸に,X[i] を縦軸に plot しますが,なんか表記のしかたが難しくて覚えにくいです(個人の感想です)。

In [16]:
plot X using (T[$1]):2 pt 7  notitle
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 7 8 9 10 gnuplot_plot_1

分割数 N を変えて結果をみます。分割数を変えても変わらない桁数までをもって,精度の良い数値解とします。

In [17]:
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]
}
Out[17]:
200 分割時の数値    10.0  0.99959156737865
400 分割時の数値    10.0000000000001  0.999591567508863
600 分割時の数値    10.0000000000001  0.999591567515716
800 分割時の数値    9.99999999999997  0.999591567516863

参考:解析解との比較

数値解析の本来の目的は,解析的に解けない問題を何とかして解きたい,解いたときの精度もちゃんと把握したい,ということです。

なので,上記の数値解がどの程度の精度で精確であるかを調べるために解析解と比較するのは本来は邪道なのですが,たまたま,上記の微分方程式は解析的に解けて,以下のようになります。

$$x(t) = \frac{\exp(t)}{9 + \exp(t)}$$

2階常微分方程式

次のような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) $$

を,

  • $a , b $ にいろいろな値を入れて
  • 初期条件 $t_0 = 0$ で $x_0 = x(t_0) = 3, v_0 = v(t_0) = 0$,
  • $t = t_0 = 0$ から $t_1 = 20$ までを $ N $ 分割

して解きます。

$N$ の値をいくつか変えてみて,精度を確認します。

In [18]:
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
}
Out[18]:
0.383972022819058  N = 100
0.38396728833795  N = 200
0.383966889383019  N = 400
0.383966861347239  N = 800

$N = 400$ で小数点以下7桁くらいの精度が出ているようですので,$N=400$ として続けます。

まずは,$a = 0, b = 0$,すなわち単振動の場合の方程式を解き,結果をファイルに保存します。

In [1]:
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]
}

計算結果をファイルに書き出します。

In [2]:
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
Out[2]:
# a = 0  b = 0
# t       x               v
0.05  2.99625078125  -0.1499375
0.1  2.985012496745  -0.299500234342448
0.15  2.96631323634216  -0.448314374121273
0.2  2.94019973845087  -0.596007961526841
0.25  2.90673727321101  -0.742211839546665
0.3  2.86600947935165  -0.88656057466559
0.35  2.81811815513779  -1.02869337025917
0.4  2.76318300392773  -1.16825496839724

ファイル rk2-a00-b00.dat のデータを plot します。

In [3]:
set grid
set zeroaxis
plot "rk2-a00-b00.dat" w l title "a=0, b=0"
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 -3 -2 -1 0 1 2 3 0 2 4 6 8 10 12 14 16 18 20 22 a=0, b=0 a=0, b=0

次の例として,$a = 0.5, b = 0.2$ の場合を解き,結果をファイルに数値データとして保存します。

In [4]:
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]
}
In [5]:
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つの数機計算の結果を重ねてグラフ化します。

In [6]:
replot "rk2-a05-b02.dat" w l title "a=0.5, b=0.2"
Gnuplot Produced by GNUPLOT 5.4 patchlevel 0 -3 -2 -1 0 1 2 3 0 2 4 6 8 10 12 14 16 18 20 22 a=0, b=0 a=0, b=0 a=0.5, b=0.2 a=0.5, b=0.2