gnuplotで微分方程式を解く

初期条件型の常微分方程式はgnuplot4.3以降の機能「累次評価演算子」(,)と擬似ファイル"+"を用いることで解くことができます。

Gauss法による1次微分方程式の数値解

まずは簡単に、Gauss法によって数値解を求めるやり方を説明します。Gauss法は直感的でわかりやすいのですが、誤差が大きいのが欠点です。

以下のような常微分方程式を解く、すなわちtの関数x(t)の関数形を求める、ことを考えます。

ここでF(x,t)は既知のある関数とします。Gauss法では、この方程式を単純に以下のように離散化します。

この離散化を用いて、時刻tでのxの値から、次のステップt+δtのxの値を近似的に得ることができます。Gauss法についての詳細は常微分方程式についての書籍やウェブページを参照してください。

Gauss法で微分方程式を初期条件x(t_ini)=x_iniの下で解くためには以下のようなスクリプトを用います。

set samples 1000
set xrange [t_ini:t_fin]
plot x=x_ini,t=t_ini, "+" using 1:(dt=$1-t,dx=F(x,t)*dt,t=$1,x=x+dx,x)

plotコマンドの中身を順に説明していきましょう。まず、x=x_iniとt=t_iniは変数の初期化です。 次に、擬似ファイル"+"を用いています。擬似ファイルは「set xrangeで指定した範囲にset samplesで指定した数の等間隔のデータが1列目に並んでいるファイル」として扱うことができます。最後にusingの中で累次評価演算子(,)を使って微分方程式を解いています。この演算子は,で区切られた代入を順に行っていき、最後の式の値を返すというものです。上記例では4つの代入演算を行い、最後にxの値を演算子の計算結果として返しています。最初の代入演算は時間ステップdtの値を取得します。次の演算ではGauss法を用いてxの増分dxを求めます。そしてtとxの値を更新し、最後にxの値を累次評価演算子全体の値として返しています。このplotコマンドを使うことで、x軸にtの値が、y軸に微分方程式の解x(t)の値がプロットされます。

次に具体例を示しましょう。以下のような方程式を解くことを考えます。

ただし初期条件は以下のとおりです。

これは充電されていないコンデンサーをt=0.0で電源のスイッチを入れて充電するときの方程式に相当します。この方程式の解析的な解は

となります。この方程式をGauss法で解いた結果と、解析的な解をプロットするスクリプトとその出力例は以下のとおりです。解析解と一致する解が得られていることが判ります。

x_ini=0.0    # xの初期値
t_ini=-2.0   # tの初期値
t_fin=20.0   # tの最終値

F(x,t)=(x0(t)-x)/tau   # 一次常微分方程式の右辺 dx/dt=F(x,t)
tau=1.0                # 
x0(t)=(t < 0.0 ? 0.0 : 1.0)

exact(t)=(t<0.0 ? 0.0 : 1.0-exp(-t/tau))  # 解析解

set xrange [t_ini:t_fin]
set yrange [-0.1:1.1]
set xlabel "time"
set ylabel "x(t)"
set key bottom inside samplen 0.8
set samples 1000
plot x=x_ini,t=t_ini, "+" u 1:(dt=$1-t,dx=F(x,t)*dt,t=$1,x=x+dx,x)\
     w p title "Gauss method",\
     exact(x) w l lw 2 title "exact solution"

Runge-Kutta法を用いた1次微分方程式の数値解

Runge-Kutta法は、Gauss法よりも数段誤差が少ない優れた解法として知られています。よく知られているものは4次のRunge-Kutta法と呼ばれているもので、時刻tでのx(t)の値から以下のように次の時刻t+δtでのx(t+δt)を求めます。

この解法による解も、以下のようなスクリプトを用いて得ることができます。

N=1000  # 分割数
dt=0.01  # 時間ステップの幅
set samples N            # 擬似ファイルのデータ数
set xrange [0:dt*(N-1)]  # 擬似ファイルのレンジ
plot xold=x_ini,\
    "+" using ($1+dt):\
    (z1=F(xold,$1), \
     z2=F(xold+z1*dt*0.5,$1+dt*0.5), \
     z3=F(xold+z2*dt*0.5,$1+dt*0.5), \
     z4=F(xold+z3*dt,$1+dt), \
     xnew=xold+(dt/6.0)*(z1+2.0*z2+2.0*z3+z4), \
     xold=xnew, \
     xnew)

実行例としては、カオスで有名なローレンツアトラクタを示します。ローレンツアトラクタの方程式は

と表され、Gauss法では上手く解けないらしいですが、4次のRunge-Kutta法ならば解くことができます。なお、この例では3次元の常微分方程式を解いているのですが、1次元のRunge-Kutta法から簡単に拡張することができます。

# Solve a set of differential equations
#  dx/dt = Fx(x,y,z,t)
#  dy/dt = Fy(x,y,z,t)
#  dz/dt = Fz(x,y,z,t)

# ローレンツアトラクタの微分方程式
Fx(x,y,z,t)=-p*x+p*y
Fy(x,y,z,t)=-x*z+r*x-y
Fz(x,y,z,t)=x*y-b*z

# ローレンツアトラクタのパラメータ
p=10.0
r=28.0
b=8.0/3.0

# ステップ数と時間ステップ幅
N=10000
dt = 0.002

# 粒子初期位置
x0=10.0
y0=20.0
z0=30.0

# Runge-Kuttaステップの変数
Dx1=0.0
Dx2=0.0
Dx3=0.0
Dx4=0.0

Dy1=0.0
Dy2=0.0
Dy3=0.0
Dy4=0.0

Dz1=0.0
Dz2=0.0
Dz3=0.0
Dz4=0.0

set samples N
set xrange [0:dt*(N-1)]

outputfile="Lorentz-Attractor.dat"
set table outputfile      # いったん解を別のファイルに書き出す

plot xold=x0,yold=y0,zold=z0,\
     "+" using ($1+dt):\
     (Dx1= Fx(xold,yold,zold,$1), \
     Dy1 = Fy(xold,yold,zold,$1), \
     Dz1 = Fz(xold,yold,zold,$1), \
     Dx2 = Fx(xold+Dx1*dt*0.5, yold+Dy1*dt*0.5, zold+Dz1*dt*0.5, $1+dt*0.5), \
     Dy2 = Fy(xold+Dx1*dt*0.5, yold+Dy1*dt*0.5, zold+Dz1*dt*0.5, $1+dt*0.5), \
     Dz2 = Fz(xold+Dx1*dt*0.5, yold+Dy1*dt*0.5, zold+Dz1*dt*0.5, $1+dt*0.5), \
     Dx3 = Fx(xold+Dx2*dt*0.5, yold+Dy2*dt*0.5, zold+Dz2*dt*0.5, $1+dt*0.5), \
     Dy3 = Fy(xold+Dx2*dt*0.5, yold+Dy2*dt*0.5, zold+Dz2*dt*0.5, $1+dt*0.5), \
     Dz3 = Fz(xold+Dx2*dt*0.5, yold+Dy2*dt*0.5, zold+Dz2*dt*0.5, $1+dt*0.5), \
     Dx4 = Fx(xold+Dx3*dt, yold+Dy3*dt, zold+Dz3*dt, $1+dt), \
     Dy4 = Fy(xold+Dx3*dt, yold+Dy3*dt, zold+Dz3*dt, $1+dt), \
     Dz4 = Fz(xold+Dx3*dt, yold+Dy3*dt, zold+Dz3*dt, $1+dt), \
     xnew=xold+(dt/6.0)*(Dx1+2.0*Dx2+2.0*Dx3+Dx4), \
     ynew=yold+(dt/6.0)*(Dy1+2.0*Dy2+2.0*Dy3+Dy4), \
     znew=zold+(dt/6.0)*(Dz1+2.0*Dz2+2.0*Dz3+Dz4), \
     xold=xnew, \
     yold=ynew, \
     zold=znew, \
     xnew):(ynew):(znew):(0):(0) w xyerrorbars

unset table

set autoscale x
set autoscale y
set autoscale z

set view 79,195
splot outputfile using 2:3:4 w l title "Lorentz Attractor"

このスクリプトでは、データの列数が多いため、set tableコマンドを使って一旦データを別ファイルに書き出しています。上記例の実行結果は以下のようになります。

このように少し複雑ではありますが、gnuplotのスクリプトだけでローレンツアトラクタをプロットすることができました。初期値を少しずらすだけで軌跡が変わることを確認してみてください。

Runge-Kutta法を用いた2次微分方程式の数値解

n元の2次微分方程式は2n元の1次微分方程式と読み替えて解くことができます。例えば、xについての2次微分方程式はxとv=dx/dtについての連立1次微分方程式と考えることができます。

以下に例を示します。この例では以下のような微分方程式を解いています。

この方程式は距離の2乗に比例する中心力の下で運動する粒子に少しだけ速度に比例する粘性項が加わった場合の運動方程式です。

# Solve a set of second-order differential equations
#  d2x/dt2 = Fx(x,y,dx/dt,dy/dt,t)
#  d2y/dt2 = Fy(x,y,dx/dt,dy/dt,t)
#  

# 微分方程式の右辺
Fx(x,y,vx,vy,t)=(-1.0/(x**2+y**2))*cos(atan2(y,x))-0.01*vx
Fy(x,y,vx,vy,t)=(-1.0/(x**2+y**2))*sin(atan2(y,x))-0.01*vy

# 変更不要
Fvx(x,y,vx,vy,t)=vx
Fvy(x,y,vx,vy,t)=vy

# Runge-Kuttaステップの時間幅とステップ数
dt = 0.01
N  = 5000

# x,yの初期値
x0=1.0
y0=0.0

# vx=dx/dt, vy=dy/dtの初期値
vx0=0.0
vy0=1.1

# Runge-Kutta法の変数の初期化
Dx1=0.0; Dx2=0.0; Dx3=0.0; Dx4=0.0
Dy1=0.0; Dy2=0.0; Dy3=0.0; Dy4=0.0
Dvx1=0.0; Dvx2=0.0; Dvx3=0.0; Dvx4=0.0
Dvy1=0.0; Dvy2=0.0; Dvy3=0.0; Dvy4=0.0

outputfile="differential-eq-output.dat"

set samples N
set xrange [0:dt*(N-1)]

# いったん別ファイルにデータを書き出す
set table outputfile
plot xold=x0,yold=y0,vxold=vx0,vyold=vy0,\
     "+" using ($1+dt):\
     (Dx1=Fvx(xold,yold,vxold,vyold,$1), \
     Dy1 =Fvy(xold,yold,vxold,vyold,$1), \
     Dvx1=Fx(xold,yold,vxold,vyold,$1), \
     Dvy1=Fy(xold,yold,vxold,vyold,$1), \
     Dx2 =Fvx(xold+Dx1*dt*0.5, yold+Dy1*dt*0.5, vxold+Dvx1*dt*0.5, vyold+Dvy1*dt*0.5, $1+dt*0.5), \
     Dy2 =Fvy(xold+Dx1*dt*0.5, yold+Dy1*dt*0.5, vxold+Dvx1*dt*0.5, vyold+Dvy1*dt*0.5, $1+dt*0.5), \
     Dvx2 =Fx(xold+Dx1*dt*0.5, yold+Dy1*dt*0.5, vxold+Dvx1*dt*0.5, vyold+Dvy1*dt*0.5, $1+dt*0.5), \
     Dvy2 =Fy(xold+Dx1*dt*0.5, yold+Dy1*dt*0.5, vxold+Dvx1*dt*0.5, vyold+Dvy1*dt*0.5, $1+dt*0.5), \
     Dx3 = Fvx(xold+Dx2*dt*0.5, yold+Dy2*dt*0.5, vxold+Dvx2*dt*0.5, vyold+Dvy2*dt*0.5, $1+dt*0.5), \
     Dy3 = Fvy(xold+Dx2*dt*0.5, yold+Dy2*dt*0.5, vxold+Dvx2*dt*0.5, vyold+Dvy2*dt*0.5, $1+dt*0.5), \
     Dvx3 = Fx(xold+Dx2*dt*0.5, yold+Dy2*dt*0.5, vxold+Dvx2*dt*0.5, vyold+Dvy2*dt*0.5, $1+dt*0.5), \
     Dvy3 = Fy(xold+Dx2*dt*0.5, yold+Dy2*dt*0.5, vxold+Dvx2*dt*0.5, vyold+Dvy2*dt*0.5, $1+dt*0.5), \
     Dx4 = Fvx(xold+Dx3*dt, yold+Dy3*dt, vxold+Dvx3*dt, vyold+Dvy3*dt, $1+dt), \
     Dy4 = Fvy(xold+Dx3*dt, yold+Dy3*dt, vxold+Dvx3*dt, vyold+Dvy3*dt, $1+dt), \
     Dvx4 = Fx(xold+Dx3*dt, yold+Dy3*dt, vxold+Dvx3*dt, vyold+Dvy3*dt, $1+dt), \
     Dvy4 = Fy(xold+Dx3*dt, yold+Dy3*dt, vxold+Dvx3*dt, vyold+Dvy3*dt, $1+dt), \
     xnew=xold+(dt/6.0)*(Dx1+2.0*Dx2+2.0*Dx3+Dx4), \
     ynew=yold+(dt/6.0)*(Dy1+2.0*Dy2+2.0*Dy3+Dy4), \
     vxnew=vxold+(dt/6.0)*(Dvx1+2.0*Dvx2+2.0*Dvx3+Dvx4), \
     vynew=vyold+(dt/6.0)*(Dvy1+2.0*Dvy2+2.0*Dvy3+Dvy4), \
     xold=xnew, \
     yold=ynew, \
     vxold=vxnew, \
     vyold=vynew, \
     xnew):(ynew):(vxnew):(vynew):(0) with xyerrorbars

unset table

set xrange [-1.5:1.5]
set yrange [-1.5:1.5]
set size ratio 1

v0=8.0
plot outputfile u 2:3 w l title "motion",\
     outputfile u 2:3:($4/v0):($5/v0) every 100 w vector lw 2 title "velocity"

このスクリプトの実行結果は以下の様になります。粒子の軌跡に加え、100点おきに粒子の速度ベクトルも表示させてみました。粘性項によって徐々に粒子の軌跡の半径が減少していくことがわかります。

おまけ

以下のスクリプトでローレンツアトラクタのアニメーションを見られます。面白いのでやってみてください。

# Solve a differential equation
#  dx/dt = Fx(x,y,z,t)
#  dy/dt = Fy(x,y,z,t)
#  dz/dt = Fz(x,y,z,t)

Fx(x,y,z,t)=-p*x+p*y
Fy(x,y,z,t)=-x*z+r*x-y
Fz(x,y,z,t)=x*y-b*z

p=10.0
r=28.0
b=8.0/3.0

dt = 0.002
N=10000

x0=10.0
y0=20.1
z0=30.0

Dx1=0.0;Dx2=0.0;Dx3=0.0;Dx4=0.0
Dy1=0.0;Dy2=0.0;Dy3=0.0;Dy4=0.0
Dz1=0.0;Dz2=0.0;Dz3=0.0;Dz4=0.0

outputfile=sprintf("Lorentz-Attractor_p%.2f_r%.2f_b%.2f.dat", p, r, b)

if (exist("i")==0 || i<0) i=0
if (i==0) set samples N
if (i==0) set xrange [0:dt*(N-1)]
if (i==0) set table outputfile

if (i==0) plot xold=x0,yold=y0,zold=z0,\
     "+" using ($1+dt):\
     (Dx1= Fx(xold,yold,zold,$1), \
     Dy1 = Fy(xold,yold,zold,$1), \
     Dz1 = Fz(xold,yold,zold,$1), \
     Dx2 = Fx(xold+Dx1*dt*0.5, yold+Dy1*dt*0.5, zold+Dz1*dt*0.5, $1+dt*0.5), \
     Dy2 = Fy(xold+Dx1*dt*0.5, yold+Dy1*dt*0.5, zold+Dz1*dt*0.5, $1+dt*0.5), \
     Dz2 = Fz(xold+Dx1*dt*0.5, yold+Dy1*dt*0.5, zold+Dz1*dt*0.5, $1+dt*0.5), \
     Dx3 = Fx(xold+Dx2*dt*0.5, yold+Dy2*dt*0.5, zold+Dz2*dt*0.5, $1+dt*0.5), \
     Dy3 = Fy(xold+Dx2*dt*0.5, yold+Dy2*dt*0.5, zold+Dz2*dt*0.5, $1+dt*0.5), \
     Dz3 = Fz(xold+Dx2*dt*0.5, yold+Dy2*dt*0.5, zold+Dz2*dt*0.5, $1+dt*0.5), \
     Dx4 = Fx(xold+Dx3*dt, yold+Dy3*dt, zold+Dz3*dt, $1+dt), \
     Dy4 = Fy(xold+Dx3*dt, yold+Dy3*dt, zold+Dz3*dt, $1+dt), \
     Dz4 = Fz(xold+Dx3*dt, yold+Dy3*dt, zold+Dz3*dt, $1+dt), \
     xnew=xold+(dt/6.0)*(Dx1+2.0*Dx2+2.0*Dx3+Dx4), \
     ynew=yold+(dt/6.0)*(Dy1+2.0*Dy2+2.0*Dy3+Dy4), \
     znew=zold+(dt/6.0)*(Dz1+2.0*Dz2+2.0*Dz3+Dz4), \
     xold=xnew, \
     yold=ynew, \
     zold=znew, \
     xnew):(ynew):(znew):(0):(0) with xyerrorbars


if (i==0) unset table

if (i==0) set autoscale xy
if (i==0) splot outputfile u 2:3:4 w l lc rgb "white"

if (i==0) set xrange [GPVAL_X_MIN:GPVAL_X_MAX]
if (i==0) set yrange [GPVAL_Y_MIN:GPVAL_Y_MAX]
if (i==0) set zrange [GPVAL_Z_MIN:GPVAL_Z_MAX]


if (i==0) set view 79,195
splot outputfile u 2:3:4 every ::::i w l title "", \
     "" u 2:3:4 every ::i::i w p pt 7 ps 2 title ""

i=i+50

if (i<N) reread

i=-1