processingの備忘録　-畑政義写像- - プログラミングの備忘録

こんにちは。
今回は「畑政義写像」を描画してみます。

そもそもは二重振り子を実装しようと思っていろいろ調べていたのですが、そのときにこんなサイトを見つけました。
（その後、二重振り子は無事実装できました ( processingの備忘録　-二重振り子- )。）

shadowacademy.web.fc2.com

その中の「非線形物理学（カオス理論）」のところに「グモウスキー・ミラの写像」というものがあり、どういうものなのかと思って調べてみたら一目惚れです。

かっこいい（小並感） pic.twitter.com/pthLxIJAkI
— taq (@taq2777) 2022年9月25日

この関連で、写像ないしはストレンジアトラクタを漁るようになり、そこで見つけたのが「畑政義写像」でした。
（少し前の話だったので、どこでこの写像を見つけたのかはわかりませんでした…）

畑政義写像とは
プログラム化
まとめ
参考

畑政義写像とは

畑政義写像は、以下のような式で表されるものです。
（写像ってなんすか…？って感じかもしれませんが、語の意味はそんなに関係ないので気になる方は調べてみてください。）

$$ \begin{align} f_{1}(z) &= a z + b \bar{z} \\ f_{2}(z) &= c (z - 1)+d (\bar{z} - 1) + 1 \end{align} $$

ここで、$a, b, c, d$ は複素数の定数、$z$ は複素数 $x + y i$、$\bar{z}$ は $z$ の共役複素数 $x - y i$ を表します。

実際に描画するときは、初期値 $z_{0}$ に対して

$$ \begin{align} z_{n+1} &= f_{1}(z_{n}) \\ z_{n+1} &= f_{2}(z_{n}) \end{align} $$

というように計算を行い、得られた $z_{n+1}$ をプロットしていきます。
そのため描画すべき点は $2^{n}$ になります。

プログラム化

数学的には複素数を扱うので、描画は複素平面に行います。
ですが、わざわざ複素数で考えなくても、実部と虚部で分けて考えれば良いです。

各複素数を $a_{ \mathrm{R} } + a_{ \mathrm{I} } i$ の形で表し、先程示した式に代入して変形していくと、以下のような漸化式が得られます。
（上のものが $z_{n+1} = f_{1}(z_{n})$、下のものが $z_{n+1} = f_{2}(z_{n})$ に対応しています。）

$$ \begin{align} x_{n+1} &= a_{ \mathrm{R} } x_{n} - a_{ \mathrm{I} } y_{n} + b_{ \mathrm{R} } x_{n} + b_{ \mathrm{I} } y_{n} \\ y_{n+1} &= a_{ \mathrm{R} } y_{n} + a_{ \mathrm{I} } x_{n} - b_{ \mathrm{R} } y_{n} + b_{ \mathrm{I} } x_{n} \end{align} $$

$$ \begin{align} x_{n+1} &= c_{ \mathrm{R} } (x_{n}-1) - c_{ \mathrm{I} } y_{n} + d_{ \mathrm{R} } (x_{n}-1) + b_{ \mathrm{I} } y_{n} + 1 \\ y_{n+1} &= c_{ \mathrm{R} } y_{n} + c_{ \mathrm{I} } (x_{n}-1) - d_{ \mathrm{R} } y_{n} + d_{ \mathrm{I} } (x_{n}-1) \end{align} $$

これで実装すべき式が得られました。

あとはこれを組み込むだけです。

ArrayList<PVector> pos = new ArrayList<PVector>();
float aR, aI, bR, bI, cR, cI, dR, dI;
float scale;
int n;

void setup(){
  size(700, 700);
  
  scale = 400; // スケール
  n = 20; // 繰り返し回数
  
  // 初期値
  PVector z0 = new PVector(0, 0);
  pos.add(z0);
  
  // パラメータ
  //aR = 0.4614; aI = 0.4614;
  //bR = 0;      bI = 0;
  //cR = 0.622;  cI = -0.196;
  //dR = 0;      dI = 0;
}

void draw(){
  background(0);
  translate(width/2, height/2);
  
  // 点を打つ座標の計算
  ArrayList<PVector> next_pos = new ArrayList<PVector>();
  for(PVector p: pos){
    next_pos.add(new PVector(aR*p.x - aI*p.y + bR*p.x + bI*p.y,
                             aR*p.y + aI*p.x - bR*p.y + bI*p.x));
    next_pos.add(new PVector(cR*(p.x-1) - cI*p.y + dR*(p.x-1) + dI*p.y + 1,
                             cR*p.y + cI*(p.x-1) - dR*p.y + dI*(p.x-1)));
  }
  pos = next_pos;
  
  // 点を打つ
  for(PVector p: pos){
    stroke(255);
    point(scale*p.x, scale*p.y);
  }
  
  n --;
  
  // 繰り返し数が上限に達したら終了
  if(n <= 0){
    noLoop();
  }
}