NyaplotでStrange Attractorを作る

## はじめに
[アトラクター](http://ja.wikipedia.org/wiki/%E3%82%A2%E3%83%88%E3%83%A9%E3%82%AF%E3%82%BF%E3%83%BC)は物理の世界から出てきたトピックで、今も活発に研究が続いているようです。その中でもストレンジアトラクターと言われるものは身近な存在で、見た目が面白いのでスクリーンセーバーなど色々なところでみるかと思います。この記事では[nyaplot](https://github.com/domitry/nyaplot)とrb-gslを使ってStrange Attractorの内いくつかを作ってみます。Nyaplotについてはこの[ブログ記事](http://domitry.hatenablog.jp/entry/2014/08/23/215630)をご覧ください。
記事中のコードは全て[こちらのノートブック](http://nbviewer.ipython.org/github/domitry/nyaplot-notebooks/blob/master/Strange_Attractors.ipynb)でご覧いただけます。

## 準備
必要なgemをインストールします。nyaplotは`gem install nyaplot`でインストールできます。必須ではありませんが`gem install nyaplot-utils`でnyaplotの補助ライブラリがインストールできます。
また、rb-gslのインストール方法については[こちら](http://qiita.com/domitry/items/e0782628699bfcadfe05)、IRubyのインストール方法については[こちら](http://qiita.com/domitry/items/aa348caf993e876a7be5)をご覧ください。
インストール後全てのライブラリをrequireします。

```rb
require 'gsl'
require 'nyaplot'
require 'nyaplot3d'
require 'nyaplot_utils'
```

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

```math
\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はこの式の一部を改変して渦が複数になるようにしたもので、その式は以下のようになります。

```math
\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モジュールのドキュメントは[こちら](http://blackwinter.github.io/rb-gsl/rdoc/odeiv_rdoc.html)。

```rb
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`に入れて内容を確認してみます。

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

IRuby上であれば以下のようなテーブルが表示されるはずです。
![chua_table.png](https://qiita-image-store.s3.amazonaws.com/0/32635/763177f3-7dc5-ce76-faa3-fa0230b6860b.png)


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

![chua_modified.png](https://qiita-image-store.s3.amazonaws.com/0/32635/99c1b267-36c5-94e0-5571-03587c652f80.png)

```rb
plot = Nyaplot::Plot.new
line = plot.add_with_df(df, :line, :x, :y)
line.stroke_width(1)
plot
```

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

```rb
plot = Nyaplot::Plot3D.new
plot.add_with_df(df, :line, :x, :y, :z)
plot
```

![chua_modified_3d.png](https://qiita-image-store.s3.amazonaws.com/0/32635/279d4221-abe7-1fc6-6d71-585b4a5f85c3.png)

以上のプロットを拡大/縮小・回転などしながらご覧になるには[こちら](http://nbviewer.ipython.org/github/domitry/nyaplot-notebooks/blob/master/Strange_Attractors.ipynb#Modified-Chua's-attractor)をご覧ください。

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

```math
\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モジュールを使って解いてみます。

```rb
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](https://qiita-image-store.s3.amazonaws.com/0/32635/801982e3-e9bb-00f6-eabd-9c156b3c484f.png)

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

```rb
zplot = Nyaplot::Plot3D.new
plot.add_with_df(df, :line, :x, :y, :z)
plot
```

![tamari_3d.png](https://qiita-image-store.s3.amazonaws.com/0/32635/df5c8049-0f9d-5ab0-ed3f-c60a140350f5.png)

## 画像のエクスポート

[nyaplot-utils](https://github.com/domitry/nyaplot-utils)を使うことで作ったプロットをSVGファイルに保存することができます。

`Nyaplot::Plot#export_svg`で軸も含めた画像全体を、`Nyaplot::Plot#export_contents`でプロットの内容だけを保存することができます。
下図は`export_contents`で出力したSVGファイルをpngに変換したものです。

```rb
plot.export_contents 'context.svg'
```

![context.png](https://qiita-image-store.s3.amazonaws.com/0/32635/637557c3-4cea-838e-b018-c13482c09cc0.png)

## 雑感

色々初期値をいじってみたのですがChua attractorがあまりきれいにならなかったので上手い初期値を教えていただけると嬉しいです。
例を作るため見栄えの良いattractorがないかと探して作ったものなので正確さは期待しないようお願いします。
見栄えとプロットのしやすさを重視するならばNyaplotは使い勝手がよく、現時点において最適解なのではないかと思います。

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

- [Nyaplot notebooks](https://github.com/domitry/nyaplot-notebooks)
- [Nyaplot Tutorials and Case studies](http://nbviewer.ipython.org/github/domitry/nyaplot/blob/master/examples/notebook/Index.ipynb)