ローレンツアトラクタ

微分方程式の項で説明したローレンツアトラクタの計算方法と乱数発生関数rand(0)を利用して、星の海に浮かぶローレンツアトラクタを描きました。
#version 4.3Sep08gp
# 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=30000

# 粒子初期位置
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"

if (exist("TERM")==0 || TERM eq "WIN") set table outputfile
if (exist("TERM")==0 || TERM eq "WIN") 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

if (exist("TERM")==0 || TERM eq "WIN") unset table


set tmargin 1
set bmargin 1
set rmargin 1
set lmargin 1

set border 0
unset xtics 
unset ytics
unset ztics
unset colorbox
unset key

set multiplot

set palette model HSV functions gray,1,1 # HSV color space

set view 79,195


set object 1 rect from screen 0, 0 to screen 1, 1 back
set object 1 rect fc rgb "black" fillstyle solid 1.0

set style fill transparent solid 0.4

set samples 100
set xrange [0:1]
set yrange [0:1]
set cbrange [0:1]
plot "+" using (rand(0)):(rand(0)):(rand(0)*0.01):(rand(0)) w circles lc palette

unset object

set autoscale x
set autoscale y
set autoscale z
set autoscale cb

splot outputfile using 2:3:4:1 w l lc palette title "Lorentz Attractor"

unset multiplot

if (exist("TERM")==0 || TERM eq "WIN") pause -1;\
     TERM = "PNG";\
     set term png interlace size 500,400 ;\
     set out "Lorentz-Attractor.png";\
     reread

set term win
set out