Watanabe-Strogatz変換の可視化

## Watanabe-Strogatz変換とは
Watanabe-Strogatz変換は, 平均場を通じて相互作用する$N$個の振動子からなる大域結合ネットワークの$N$次元ダイナミクスを, 3次元の運動に縮約する方法である.

## N次元大域結合振動子ネットワーク
次の方程式で表される$N$振動子の大域結合ネットワークを考える.

$$\dot{\theta_i} = \omega + \frac{\epsilon}{N} \sum_{j=1}^{N}\cos(\theta_i - \theta_j)  \tag{1}$$
$\theta_i (i = 1,2, \dots,N)$は$i$番目の振動子の位相, $\omega$は共通の自然振動数である. 振動子同士は結合強度$\epsilon$で相互作用している.

式$(1)$を$\theta_i^{\prime}=\theta_i - \omega t$という変換によって速度$\omega$で回転する座標系に移行する. さらに, 時間を$t^{\prime}=\epsilon t$とリスケールすると, $\omega$も$\epsilon$も含まない単純な式にできる.
$$\dot{\theta_i} = \frac{1}{N} \sum_{j=1}^{N}\cos(\theta_i - \theta_j) \tag{2}$$

この式$(2)$を出発点とする.

次の図とアニメーションは$N=20$としたときのダイナミクスの一例である. この系が複雑な挙動を示すことがわかる.

![theta_ts.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/65132/69e67edd-9631-533b-5968-8ae3f20bd0f8.jpeg)

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/65132/ff883d84-01e3-030c-89ea-eba85f21e0a2.gif">

## 変換の概要
Watanabe-Strogatz変換は以下3ステップの座標変換を行なう.
1. 平均位相$\Theta(t)$とともに回転する座標系に変換（$\theta \rightarrow \xi$）
2. 円周上に偏って分布する振動子を均一に分布させる座標系に変換（$\xi \rightarrow \eta$）
3. 全ての振動子が静止するような座標系に変換（$\eta \rightarrow \psi$）

全ての座標変換をまとめると次のような式になる.

$$\tan\frac{\theta_i(t) - \Theta(t)}{2}=\sqrt\frac{1 + \gamma(t)}{1 - \gamma(t)}\tan\frac{\psi_i - \Psi(t)}{2} \tag{3}$$

式$(3)$にあらわれる3つの変数$\gamma(t)$, $\Theta(t)$, $\Psi(t)$の時間発展と$N$個の定数$\lbrace\psi_i\rbrace$によって, 元の$N$変数$\lbrace \theta_i(t) \rbrace$のダイナミクスを完全に再現することができる.

3つの変数の時間発展は次の微分方程式で記述される.

```math
\begin{align}
\dot{\gamma} &= \frac{(1-\gamma^2)^{\frac{3}{2}}}{N}\sum_{j}\frac{\sin(\psi_j-\Psi)}{1-\gamma\cos(\psi_j-\Psi)} \tag{4a}\\
\gamma\dot{\Psi} &= \frac{(1-\gamma^2)^{\frac{1}{2}}}{N}\sum_{j}\frac{\gamma-\cos(\psi_j-\Psi)}{1-\gamma\cos(\psi_j-\Psi)} \tag{4b}\\
\gamma\dot{\Theta} &= \frac{1}{N}\sum_{j}\frac{\gamma-\cos(\psi_j-\Psi)}{1-\gamma\cos(\psi_j-\Psi)} \tag{4c}
\end{align}
```
なお, 今回の例では式$(4)$のいずれの式の右辺にも$\Theta(t)$が含まれないことからわかるように, 3変数ではなく実質的に$\gamma$と$\Psi$の2変数のダイナミクスとなっており, $\Theta$はそれらに従属する.

以下に, それぞれの座標変換が何をしているのかをアニメーションで確認していく.

### Step1: 平均位相とともに回転する座標系に変換

まず, 元の運動$\lbrace\theta_i(t)\rbrace$を,平均位相$\Theta(t)$との相対運動を表す$\lbrace \xi_i(t) \rbrace$に座標変換する.

$$\xi_i(t) = \theta_i(t) - \Theta(t) \tag{4}$$

下のアニメーションでは, $\lbrace\theta_i(t)\rbrace$と$\Theta(t)$との関係がわかりやすいように振動子の集団を幅の狭い正規分布で初期化している. $\Theta(t)$は振動子の集団の中心とおおよそ$\pi$の位相差を保って回転し, したがって座標系$\lbrace\xi_i(t)\rbrace$でみると集団の位置がほぼ$\pi$のあたりに固定されることがわかる. ただし, 集団を構成する個々の振動子は静止しているわけではない. 集団の中心近くでは遅く, 集団から離れると速く, 一定方向に運動をしている.

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/65132/27857096-902c-38dc-9878-c522ecb0ce3b.gif">

### Step2: 円周上に偏って分布する振動子を均一に分布させる座標系に変換

$\lbrace\xi_i(t)\rbrace$の座標系では振動子の位置に偏りがあり, また, 各々の振動子は異なる速度で回転している. 次式によって$\lbrace\eta_i(t)\rbrace$に座標変換すると, 振動子が円周上にほぼ均一に分布し, 回転速度は一定となる.

$$\tan\frac{\xi_i(t)}{2}=\sqrt\frac{1 + \gamma(t)}{1 - \gamma(t)}\tan\frac{\eta_i(t)}{2} \tag{5}$$

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/65132/40cf3112-0893-cfd0-a36f-9da6877c62ee.gif">

この変換に含まれるパラメータ$\gamma(t)$は, $0\leqq\gamma(t)<1$の範囲の値をとり, $\gamma=0$なら恒等変換で, $\gamma$の値が大きくなるほど分布を大きく歪めるような変換となる.

### Step3: 全ての振動子が静止するような座標系に変換

上で述べたように, $\lbrace\eta_i(t)\rbrace$は全ての振動子が共通の速度$-\Psi(t)$で回転する. 速度$-\Psi(t)$で回転する座標系から見ると, 振動子は静止し, $\lbrace\psi_i\rbrace$は時間に依存しない定数となる.

$$\psi_i = \eta_i(t) + \Psi(t) \tag{6}$$

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/65132/02559828-7452-6e65-47fa-5917432858b3.gif">

## 初期条件の決定
さて, 実際に縮約後の方程式$(4)$を解くためには, 3つの変数の初期条件$\gamma(0)$, $\Psi(0)$, $\Theta(0)$と, $N$個の定数$\lbrace\psi_i\rbrace$を決める必要がある.

未知変数が$N+3$個であるのに対して, 式$(4)$は$N$個の方程式しかないので, 以下3つの式によって$\lbrace\psi_i\rbrace$に対称性を課す.

```math
\begin{align}
\sum_{j}\sin\psi_j=0 \tag{7a}\\
\sum_{j}\cos\psi_j=0 \tag{7b}\\
\sum_{j}\psi_j=0 \tag{7c}
\end{align}
```
これによって, $\lbrace\theta_i(0)\rbrace$から, 初期条件$\gamma(0)$, $\Psi(0)$, $\Theta(0)$と$N$個の定数$\lbrace\psi_i\rbrace$を決められるようになる.

## 低次元のダイナミクス
初期条件と定数を決めることができたので, 式$(4)$で$\gamma(t)$, $\Psi(t)$, $\Theta(t)$の時間発展を計算することができる. 下図に3変数の時系列の一例を示す.

![reduced_ts.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/65132/8207300d-4b40-eb27-7fde-a8a2bc9f71ee.jpeg)

$\gamma(t)$がほぼ$1$に近い値になっているのは元の座標系で振動子の分布を局在させたことによる. ばらつきの大きい分布の場合は$\gamma(t)$の値は小さくなる.

## 元のダイナミクスの再構成
式$(3)$を$\theta_i(t)$について整理した次の変換式によって, 低次元のダイナミクスから元の$N$変数のダイナミクスを再構成できる.

$$\theta_i(t) = 2\arctan\bigg\\{\sqrt\frac{1 + \gamma(t)}{1 - \gamma(t)}\tan\frac{\psi_i-\Psi(t)}{2}\bigg\\} + \Theta(t)\tag{8}$$

上が元の時系列, 下が再構成した時系列である.

![theta_ts_2.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/65132/d513f871-cfb8-3fed-16de-d79ec29e9d29.jpeg)

![reconstructed_theta_ts_2.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/65132/7498add2-302f-2663-71f9-042d2aedc79c.jpeg)

完璧に再現できていることがわかる.

## 参考文献
1. [S. Watanabe and S.H. Strogatz, _Phys. Rev. Lett._ **70**, 2391 (1993)](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.70.2391)
2. [S. Watanabe and S.H. Strogatz, _Physica D_ **74**, 197 (1994)](https://www.sciencedirect.com/science/article/abs/pii/0167278994901961)