Edited at

NyaplotでStrange Attractorを作る

More than 5 years have passed since last update.


はじめに

アトラクターは物理の世界から出てきたトピックで、今も活発に研究が続いているようです。その中でもストレンジアトラクターと言われるものは身近な存在で、見た目が面白いのでスクリーンセーバーなど色々なところでみるかと思います。この記事ではnyaplotとrb-gslを使ってStrange Attractorの内いくつかを作ってみます。Nyaplotについてはこのブログ記事をご覧ください。

記事中のコードは全てこちらのノートブックでご覧いただけます。


準備

必要なgemをインストールします。nyaplotはgem install nyaplotでインストールできます。必須ではありませんがgem install nyaplot-utilsでnyaplotの補助ライブラリがインストールできます。

また、rb-gslのインストール方法についてはこちら、IRubyのインストール方法についてはこちらをご覧ください。

インストール後全てのライブラリをrequireします。

require 'gsl'

require 'nyaplot'
require 'nyaplot3d'
require 'nyaplot_utils'


Modified Chua attractor

元となったChua's attractorはコイルとコンデンサ2つ、抵抗、ダイオードから成る回路に電流を流したとき、両コンデンサの両端の電位差とコイルに流れる電流の間に成り立つ関係をプロットしたものだそうです。式は以下のようになります。

\frac{dx}{dt}=\alpha [y-x-f(x)]\\

\frac{dy}{dt}=x-y+z\\
\frac{dz}{dt}=-\beta y

(出展: http://en.wikipedia.org/wiki/Chua's_circuit)

Modified Chua attractorはこの式の一部を改変して渦が複数になるようにしたもので、その式は以下のようになります。

\frac{dx}{dt}= \alpha*(y-h)\\

\frac{dy}{dt}=x-y+z\\
\frac{dz}{dt}=-\beta*y\\
h := -b*sin(\frac{\pi*x}{2*a}+d)

(出展: http://en.wikipedia.org/wiki/Multiscroll_attractor)

いずれともパラメータは3つです。それではrb-gslに用意されているOdeivモジュールで上の微分方程式を解いてみます。Odeivモジュールのドキュメントはこちら

alpha = 10.82

beta = 14.286
a = 1.3
b = 0.1
d = 0.2
func = lambda {|x| return -b*Math.sin((Math::PI*x)/(2*a)+d)}

ode_func = Proc.new { |t, y, dydt, mu|
dydt[0] = alpha*(y[1] - func.call(y[0]))
dydt[1] = y[0] - y[1] + y[2]
dydt[2] = -beta*y[1]
}

solver = GSL::Odeiv::Solver.alloc(GSL::Odeiv::Step::RK8PD, [1e-11, 0.0], ode_func, 3)

t = 0.0 # initial time
t1 = 250.0 # end time
h = 1e-11 # initial step
y = GSL::Vector.alloc([0, 0.49899, 0.2]) # initial value
mat = []

while t < t1
t, h, status = solver.apply(t, t1, h, y)
break if status != GSL::SUCCESS
mat.push([t, y[0],y[1],y[2]])
end

この結果をNyaplot::DataFrameに入れて内容を確認してみます。

mat = mat.transpose

df = Nyaplot::DataFrame.new({t: mat[0], x: mat[1], y: mat[2], z: mat[3]})
df

IRuby上であれば以下のようなテーブルが表示されるはずです。

chua_table.png

一見したところN/Aもなく、上手く計算できているようです。それではこの結果をプロットしてみます。

chua_modified.png

plot = Nyaplot::Plot.new

line = plot.add_with_df(df, :line, :x, :y)
line.stroke_width(1)
plot

できました!次にNyaplot3Dを使って3Dで表示してみます。

plot = Nyaplot::Plot3D.new

plot.add_with_df(df, :line, :x, :y, :z)
plot

chua_modified_3d.png

以上のプロットを拡大/縮小・回転などしながらご覧になるにはこちらをご覧ください。


Tamari's Attractor

国の経済をモデル化してできたattractorのようです。式は以下のようになります。

\frac{dx}{dt} = (x - ay)cos(z) - b*ysin(z)\\

\frac{dy}{dt} = (x + cy)sin(z) + d*ycos(z)\\
\frac{dz}{dt} = e + fz + g*atan(\frac{(1 - u)}{(1 - i)*xy})

(出展: http://www.bentamari.com/attractors.html)

Chua attractorと同じ要領でOdeivモジュールを使って解いてみます。

a=1.013; b=-0.021;c=0.019;d=0.96;e=0;f=0.01;g=1;u=0.05;i=0.05

ode_func = Proc.new { |t, y, dydt, mu|
dydt[0] = (y[0]-a*y[1])*Math.cos(y[2])-b*y[1]*Math.sin(y[2])
dydt[1] = (y[0]+c*y[1])*Math.sin(y[2])+d*y[1]*Math.cos(y[2])
dydt[2] = e + f*y[2] + g*Math.atan(((1-u)*y[1]) / (1-i)*y[0])
}

solver = GSL::Odeiv::Solver.alloc(GSL::Odeiv::Step::RK8PD, [1e-10, 0.0], ode_func, 3)

t = 0.0
t1 = 800.0
h = 1e-10
y = GSL::Vector.alloc([0.9,1,1])
mat = []

while t < t1
t, h, status = solver.apply(t, t1, h, y)
break if status != GSL::SUCCESS
mat.push([y[0],y[1],y[2]])
end

mat = mat.transpose
df = Nyaplot::DataFrame.new({x: mat[0], y: mat[1].map{|v| -v}, z:mat[2]})

plot = Nyaplot::Plot.new
line = plot.add_with_df(df, :line, :y, :x)
line.stroke_width(0.3)
line.color(color.to_a[5])
plot

tamari.png

3Dでも表示させてみます。

zplot = Nyaplot::Plot3D.new

plot.add_with_df(df, :line, :x, :y, :z)
plot

tamari_3d.png


画像のエクスポート

nyaplot-utilsを使うことで作ったプロットをSVGファイルに保存することができます。

Nyaplot::Plot#export_svgで軸も含めた画像全体を、Nyaplot::Plot#export_contentsでプロットの内容だけを保存することができます。

下図はexport_contentsで出力したSVGファイルをpngに変換したものです。

plot.export_contents 'context.svg'

context.png


雑感

色々初期値をいじってみたのですがChua attractorがあまりきれいにならなかったので上手い初期値を教えていただけると嬉しいです。

例を作るため見栄えの良いattractorがないかと探して作ったものなので正確さは期待しないようお願いします。

見栄えとプロットのしやすさを重視するならばNyaplotは使い勝手がよく、現時点において最適解なのではないかと思います。

Nyaplotにはこのほかにも色々な例があります。是非ご覧ください。