はじめに
gnuplot はサイエンスでも使われる有名なプロットソフトウェアです．軽量であり気軽に使うことができるので電卓として使用している人も多いのではないでしょうか．
print sin(0.1*pi)
複雑な計算となるとプロッターである gnuplot よりも python などのプログラミング言語を使うことが多いと思いますが，提供されている機能を組み合わせることでそこそこ複雑な計算もこなすことができます．そして何より計算と可視化が同時にできるということが gnuplot を電卓として使うことの一番のメリットだと思います．ここでは gnuplot を用いた積分の方法について簡単に紹介します．
目標
今回は以下のようなプロットを作製することを目標にします．このプロットはサインカーブのうえで跳ねるボールの挙動をシミュレートしたものです．gnuplot 5.0 patchlevel 0 on Ubuntu 16.04 で作成しました．
特殊ファイル名
gnuplot のコンソールで help special-filenames と入力すると plot 時に使用できる特殊なファイル名の解説が表示されます．特殊ファイル名には以下の 4 つがあります．
-
'': 直前に使用したファイルをもう一度読み込む -
'-': プロットするデータをターミナルから入力する (標準入力) -
'+': 1 次元プロットでの x 軸のサンプル点をカラム 1 に読み込む -
'++': x, y 平面でのサンプル点をそれぞれカラム 1, 2 に読み込む
ここでは '+' を使用します．このファイル名を使用すると具体的にどのようなデータが出力されるのかの例を以下に示しました．普通に sin(x) をプロットすれば表示されるものを敢えて '+' を使用して同一のものをプロットしています．
set xrange [0:20]
set yrange [-1.5:1.5]
plot sin(x) w lp, '+' u 1:(sin($1)) w p pt 6 ps 2 lw 2
このサンプル点の数は set sample によって設定可能です．デフォルトでは 100 に設定されていると思います (show sample で確認可能)．ここで (sin($1)) が sample で設定された数だけ評価されていることを覚えておいてください．
数値積分
gnuplot に限らず多くの言語では代入文では代入された値が評価値として返ります．そして gnuplot では代入文を using の内部で使うことができます．using に与えた式はサンプリングの度に評価されるため，代入文を書くことで変数を更新していくことができます．以下にサンプルを示します．
set yrange [0:250]
set xrange [0:20]
set sample 100
A = 0
dx = 20/(100.-1)
plot 0.5*x**2, '+' u 1:(A=A+$1*dx,A) w p pt 6
using の内部では A の値にサンプリングされた値とサンプリング間隔である 20÷99 の積を加えています．これは近似的に以下の積分操作に相当しています．
\sum_{n=1}^{100} x_n \frac{20}{99}
~\sim~ \int_0^x x' \,\mathrm{d}x'
~=~ \frac{1}{2}x^2
つまり gnuplot で積分ができるということです．
カンマ演算子
ところで上記の結果はどの程度正しいのでしょうか．解析的な結果である $\frac{1}{2}x^2$ との違いが気になります．A と 0.5*x**2 の差をプロットすることはできるでしょうか．
カンマ演算子は gnuplot では二項演算子のひとつ (serial evaluation) として定義されています．これは (A,B) と表記することで A を評価した後に B を評価して B の値を返すという働きをします．
A = 0
print (A = A+3, A) # 3 が出力される
そしてこの演算子は using の中でも使うことができます．この演算子を使うことで目的を達成することができます．具体的には以下の通りです．
set yrange [0.0:3.0]
set xrange [0:20]
set sample 100
A = 0
dx = 20/(100.-1)
plot '+' u 1:(A=A+$1*dx,A-0.5*x**2) w p pt 6
無事に解析解との差をプロットすることができました．なお，このカンマ演算子はいくつでも連結することができます．そのため using の内部で多少複雑な計算を実行することも可能です．
計算例
ここでは以下の運動方程式 (強制振動) に従う粒子の運動を考えます．変数の上のドットは時間微分を意味しています．
\begin{align}
\dot{y} &= v, \\
\dot{v} &= - y + 0.3\cos(0.5t)
\end{align}
初期条件として $y(0)=0$, $v(0) = 0$ としたときの解は以下のようになります．
y = 0.4\left(-\cos(t)+\cos(0.5t)\right)
とりあえず解いてみます．単純にオイラー法で計算したところ発散したので中央差分を用いて計算しています．刻み幅を十分に細かくとれば計算ができるはずです．
set yrange [-1.0:1.0]
set xrange [0:50]
set sample 500
dt = 50/(500.-1)
dydt(y,v,t) = v
dvdt(y,v,t) = -y + 0.3*cos(0.5*t)
Y = 0
V = 0
plot \
0.4*(-cos(x)+cos(x/2.0)) t 'Analytical Solution', \
'+' u ($1+dt):(\
tmpv=Vp+2*dt*dvdt(Y,V,$1),\
tmpy=Yp+2*dt*dydt(Y,V,$1),\
Vp=V,Yp=Y,V=tmpv,Y=tmpy\
) w lp pt 6 t 'Calculation'
以下のように関数を定義すれば Runge=Kutta 法で計算することもできます．ただし以下の式は陽解法なので不安定であり刻み幅を細かくすると簡単に発散します．
set yrange [-1.0:1.0]
set xrange [0:50]
set sample 10000
dt = 50/(10000-1.)
dydt(y,v,t) = v
dvdt(y,v,t) = -y + 0.3*cos(0.5*t)
rk4_dydt(y,v,t) = (\
k1=dydt(y,v,t), \
k2=dydt(y,v+dt*k1/2,t+dt/2), \
k3=dydt(y,v+dt*k2/2,t+dt/2), \
k4=dydt(y,v+dt*k3,t+dt), \
(k1+2*k2+2*k3+k4)/6.0)
rk4_dvdt(y,v,t) = (\
k1=dvdt(y,v,t), \
k2=dvdt(y+dt*k1/2,v,t+dt/2), \
k3=dvdt(y+dt*k2/2,v,t+dt/2), \
k4=dvdt(y+dt*k3,v,t+dt), \
(k1+2*k2+2*k3+k4)/6.0)
Y = Yp = 0.
V = Vp = 0.
plot \
0.4*(-cos(x)+cos(x/2.0)) t 'Analytical Solution', \
'+' u ($1+dt):(\
tmpv=V+dt*rk4_dvdt(Y,V,$1),\
tmpy=Y+dt*rk4_dydt(Y,V,$1),\
V=tmpv,Y=tmpy \
) w p pt 6 t 'Calculation'
ここからは冒頭で紹介したサンプルを作成します．重力加速度を $g$ とすると質点の従う運動方程式は以下のラグランジアンから導くことができます．
\mathcal{L}(v_x,v_y,x,y) = \frac{1}{2}(v_x^2+v_y^2) - gy
ここで $v_x,v_y,x,y$ は質点の速度及び座標を表しています．ここで空間には $y=f(x)$ で定義される壁があり，質点はこの曲線よりも下側には侵入できないとします．つまり質点の座標 $x,y$ は $y > f(x)$ という関係式を満たします．運動方程式にこのような制限を与えるためにはラグランジュの未定係数法を用います．
\mathcal{L}(v_x,v_y,x,y,\lambda) = \frac{1}{2}(v_x^2+v_y^2) - gy
+ \lambda(y-f(x)), \quad (\lambda \ge 0)
作用積分を小さくすることを考えると定性的には $y-f(x)>0$ が満たされている時には $\lambda = 0$ でありそれ以外の場合には $\lambda$ は正の値をとります．ラグランジアンから運動方程式を導出すると以下のようになります．
\begin{align}
\dot{v}_y &= -g + \lambda, \\
\dot{v}_x &= -\lambda f'(x),
\end{align}
$\lambda$ が関係している項に注目すると $\left(\lambda, -\lambda f'(x)\right)$ という力が加わっていることになります．これはちょうど $y=f(x)$ と直行する向きです．壁面からの力は常に法線の向きに働くことを考えると，$\lambda$ の関係する項は壁面からの垂直抗力を表していることがわかります．$\lambda$ は物理的には壁面から受ける力の $y$ 成分に相当します．「$y-f(x)>0$ が満たされている時には $\lambda = 0$ である」ということは物理的には壁と接触しなければ抗力が発生しないことと同義です．
このままでは $\lambda$ の値が定まりません．$\lambda$ は壁と質点がどのような相互作用をするかで決まります．ここでは完全弾性衝突をすると仮定します．すると衝突の前後で運動エネルギーが保存されなければなりません．この制約から次の方程式が得られます．
\left( v_x + \mathrm{d}t \dot{v}_x \right)^2 + \left( v_y + \mathrm{d}t \dot{v}_y \right)^2 ~=~ v_x^2 + v_y^2
この式を $\lambda$ について解くと以下の式が得られます．
\mathrm{d}t(f'(x)^2 + 1)\lambda^2
- 2(v_xf'(x) - v_y + g\mathrm{d}t)\lambda
- g^2\mathrm{d}t + 2v_yg = 0.
二次方程式なので 2 つ解がありますが符号が正のものを採用します．違う方を採用するとボールが壁にめり込んでいきます．
これで $\lambda$ の値がわかったので運動方程式を計算できるようになりました．実際に計算してみます．
set xrange [-20:20]
set yrange [-2:8]
set sample 10000
dt = 20./(10000-1.)
g = 9.8
X = pX = x0 = 0.0
Y = pY = y0 = 3.0
VX = pVX = vx0 = 3.0
VY = pVY = vy0 = 1.0
C2(x) = dt*(1+cos(x)**2)
C1(vx,vy,x) = (vx*cos(x)-vy+g*dt)
C0(vy) = g**2*dt-2*vy*g
lambda(vx,vy,x,y) = (\
(y-sin(x))>0)?(0.0):\
((C1(vx,vy,x)+sqrt(C1(vx,vy,x)**2+C2(x)*C0(vy)))/C2(x))
dvydt(vx,vy,x,y) = (-g+lambda(vx,vy,x,y))
dvxdt(vx,vy,x,y) = -lambda(vx,vy,x,y)*cos(x)
dxdt(vx,vy,x,y) = vx
dydt(vx,vy,x,y) = vy
plot \
sin(x) t "Wall" lw 2, \
'+' u (x0):(y0) w p pt 7 not,\
'+' u (\
pVX=VX+dt*dvxdt(VX,VY,X,Y),\
pVY=VY+dt*dvydt(VX,VY,X,Y), \
X=X+dt*dxdt(pVX,pVY,X,Y),\
Y=Y+dt*dydt(pVX,pVY,X,Y), \
VX=pVX,VY=pVY,X\
):(\
Y\
) w l t "Ball"
ボールの初期位置をドットで表しています．ボールの軌跡を水色の線で表示しました．とりあえずボールが跳ねているような挙動になっています．衝突判定の問題もありますし残念ながら刻み幅の変更に対してはあまり安定な解ではありませんが，それっぽいものが得られたということで満足したいと思います．
Enjoy gnuplot!
