gnuplotでアニメーション
gnuplotでルンゲクッタ法

普通、数値計算した結果をプロットする場合、C言語などで計算した結果をgnuplot などの描画用ソフトでプロットするのが一般的ですが、 手軽にちょっと確かめたい場合など、gnuplot だけで結果を確認出来れば便利です。
今回は、gnuplot をもちいて、常微分方程式の数値解を４次のルンゲクッタ法で計算し、直ちに描画するプログラムを書きます。

■参考図書：gnuplotの精義―フリーの高機能グラフ作成ツールを使いこなす [単行本]

ローレンツ方程式を解く

ローレンツ方程式についてはこちら（Wikkipedia）をご覧ください。

ローレンツアトラクタの動画

ローレンツ方程式は、時間発展がカオス的な振る舞いをすることが知られています。
各時刻ごとに jpg ファイルを書き出し、連結することで動画を作成します。
上記のことを実現するために、メインとサブルーチンの２つのgnuplotのプログラムファイルを用意します。

メインルーチン：Lorenz_Attractor.plt

gnuplot の出力の設定や方程式、パラメータの設定を行ないます。

set terminal jpeg  enhanced font "Times" 20 size 600, 480
set tics font 'Times,12'
set nokey
set ticslevel 0

set xr[-50:50]
set yr[-50:50]
set zr[0:50]

DATA = "Lorenz_Attractor.data"
set print DATA
f1(x, y, z) = -p*x + p*y
f2(x, y, z) = -x*z + r*x - y
f3(x, y, z) = x*y - b*z
p = 10.0
r = 28.0
b = 8.0/3.0
x = 1
y = 1
z = 1
dt = 0.01    #時間刻み
n_max = 1000 #静止画数
m_max = 10   #間引き数
n=0;
m=0
call "rk.plt"

計算精度を上げるために時間間隔「dt」は小さく取る必要がありますが、 全ての時刻で描画すると、時間経過を見る上で、非常に多くの枚数が必要になります。 そのために描画を間引くことで任意の枚数に仕上げることができます。 上の「m_max=10」は、数値計算10ステップごとに1回描画することを想定して、変数を定義します。

サブルーチン：rk.plt

各時刻ごとに呼び出すサブルーチン。 「reread」でサブルーチンファイルを呼び出します。

k1x = dt*f1(x, y, z)
k1y = dt*f2(x, y, z)
k1z = dt*f3(x, y, z)
k2x = dt*f1(x+k1x/2,y+k1y/2, z+k1z/2)
k2y = dt*f2(x+k1x/2,y+k1y/2, z+k1z/2)
k2z = dt*f3(x+k1x/2,y+k1y/2, z+k1z/2)
k3x = dt*f1(x+k2x/2,y+k2y/2, z+k2z/2)
k3y = dt*f2(x+k2x/2,y+k2y/2, z+k2z/2)
k3z = dt*f3(x+k2x/2,y+k2y/2, z+k2z/2)
k4x = dt*f1(x+k3x,  y+k3y,   z+k3z)
k4y = dt*f2(x+k3x,  y+k3y,   z+k3z)
k4z = dt*f3(x+k3x,  y+k3y,   z+k3z)
x=x+(k1x + 2*k2x + 2*k3x + k4x)/6
y=y+(k1y + 2*k2y + 2*k3y + k4y)/6
z=z+(k1z + 2*k2z + 2*k3z + k4z)/6
print x, y, z
m=m+1
if(m < m_max) reread
m=0

outfile(n) = sprintf("f/%d.jpg",n+10000)  #出力ファイル名
title(n) = sprintf("t = %4.2f",n*dt )  #タイトル名
unset label 
set label title(n)  font 'Times,12'  at -50 , -50 , 55 
set output outfile(n)
splot DATA with line linecolor rgb "blue"

n = n+1
if(n < n_max) reread
n=0