【OculusGo】Loranz方程式のカオスな振舞をVR♬

[先ほどDuffing振動子のカオスな振舞をVR](https://qiita.com/MuAuan/items/82dbdb673b799a615fd1)したが、今回はカオス研究の発端となったLorenz方程式のカオスな振舞をVRしようと思う。

こちらは、Duffing振動子は一次元二階常微分方程式を解いたが、今回は一階の３次元常微分方程式を解く必要がある。
ちょっと目先が異なるのでodeintで解けるものやらと思っていたが、解けたので記載しておこうと思う。

###やったこと
（１）Lorenz方程式の解法
（２）Lorenz方程式の振舞
（３）OculusGoでVRしてみる
###（１）Lorenz方程式の解と振舞
コードは以下に置いた
・[LSTM/duffing.py](https://github.com/MuAuan/LSTM/blob/master/duffing.py)
大切なことは、どちらの考え方でも解けた。
つまり、通常の運動方程式と考えて、再度微分して３次元二階常微分方程式でもいいが、普通に３次元一階常微分方程式と考えて解けた。

```py
def lorenz(var, t, p, r, b):
    """
    var = [x,y,z, x_dot,y_dot,z_dot]
    dx/dt = -px+py
    dy/dt = -xz+rx-y
    dz/dt = xy-bz
    """
    x_dot = -p*var[0]+p*var[1]
    y_dot = -var[0]*var[2]+r*var[0]-var[1]
    z_dot = var[0]*var[1]-b*var[2]
    px_dot= -p*x_dot +p*y_dot
    py_dot= -x_dot*var[2] - var[0]*z_dot +r*x_dot -y_dot
    pz_dot= x_dot*var[1] + var[0]*y_dot -b*z_dot
    return x_dot,y_dot,z_dot,px_dot,py_dot,pz_dot
```

あるいは、以下のようなものである。

```py
def lorenz(var, t, p, r, b):
    """
    var = [x,y,z]
    dx/dt = -px+py
    dy/dt = -xz+rx-y
    dz/dt = xy-bz
    """
    x_dot = -p*var[0]+p*var[1]
    y_dot = -var[0]*var[2]+r*var[0]-var[1]
    z_dot = var[0]*var[1]-b*var[2]
    return x_dot,y_dot,z_dot
```

```py
p = 10
b = 8/3
r = 28
y0 = [0,10,20]
t = np.linspace(0, 35, 3001)
sol = odeint(lorenz, y0, t, args=(p, r, b))
x, y, z = sol.T[0], sol.T[1], sol.T[2]
```
上記の二通りのやり方は、初期値を除いてほぼ同じようなものなので（二階部分は単なる検証程度の意味しかない）、実際一階の方程式が解けているにすぎなさそうである。
###（２）Lorenz方程式の振舞

```math
{\frac  {dx}{dt}}=-px+py
```

```math
{\frac  {dy}{dt}}=-xz+rx-y
```

```math
{\frac  {dz}{dt}}=xy-bz
```
結果は、以下の参考の動画像の通りだが、。。やってみると結構面白い結果になった。ここでは以下代表的と思われる初期値等のセットを掲載する。上記のコードで比較的簡単に適当な初期値で計算できるので試せると思う。
【参考】
・[ローレンツ方程式](https://ja.wikipedia.org/wiki/%E3%83%AD%E3%83%BC%E3%83%AC%E3%83%B3%E3%83%84%E6%96%B9%E7%A8%8B%E5%BC%8F)
・[File:A Trajectory Through Phase Space in a Lorenz Attractor.gif](https://commons.wikimedia.org/wiki/File:A_Trajectory_Through_Phase_Space_in_a_Lorenz_Attractor.gif)

まず、上記参考の通りの初期パラメータ
p = 10
b = 8/3
r = 28
y0 = [0,10,20]
t = np.linspace(0, 35,3001)
だと以下のとおりとなった。
plot_xy
![plot_xy.png](https://qiita-image-store.s3.amazonaws.com/0/233744/444061eb-3995-f78a-12b1-3ba7e23495a5.png)
plot_yz
![plot_yz.png](https://qiita-image-store.s3.amazonaws.com/0/233744/75678206-4438-1028-8c13-b7ad591fd5fa.png)
plot_zx
![plot_zx.png](https://qiita-image-store.s3.amazonaws.com/0/233744/195cf035-80d6-0b28-f99e-142c81a0b73a.png)
そして時間変化は
plot_xt
![plot_xt.png](https://qiita-image-store.s3.amazonaws.com/0/233744/d3f75981-67f6-981e-e30a-c1a9f77f1cce.png)
plot_yt
![plot_yt.png](https://qiita-image-store.s3.amazonaws.com/0/233744/a1453846-fd25-c150-9bb6-2f8932e917f8.png)
plot_zt
![plot_zt.png](https://qiita-image-store.s3.amazonaws.com/0/233744/a05f0ee6-5033-61b0-8448-b540ec089b5e.png)

一方、y0 = [0.1,0, 0.0]のとき
t = np.linspace(0, 20, 3001)とすると以下が得られる。
plot_xy
![plot_xy.png](https://qiita-image-store.s3.amazonaws.com/0/233744/b6d6a9af-df08-87ec-1c0a-de376b1a0ede.png)
plot_yz
![plot_yz.png](https://qiita-image-store.s3.amazonaws.com/0/233744/729e3edf-d42a-0c83-95e0-7ab991eb6a4e.png)
plot_zx
![plot_zx.png](https://qiita-image-store.s3.amazonaws.com/0/233744/85b69872-d920-2e3c-24bc-3268e035bf7d.png)
plot_xt
![plot_xt.png](https://qiita-image-store.s3.amazonaws.com/0/233744/89f267bf-b1c9-4a1c-167c-8c7130b3653a.png)
plot_yt
![plot_yt.png](https://qiita-image-store.s3.amazonaws.com/0/233744/e32db2bf-c43a-f14a-c898-246419c68e6d.png)
plot_zt
![plot_zt.png](https://qiita-image-store.s3.amazonaws.com/0/233744/1b1f9a6d-e2f0-0006-8a70-011cea471f51.png)
###（３）OculusGoでVRしてみる
さて、いよいよVRで見えるかな？？
ほぼ前回と同じなので見えるはず。。
今回は、３パラメータの授受だけど基本同じ手法でやりました。
[OculusGoでLorenz方程式のカオスな運動を見る♬](https://youtu.be/RAf4r4mEuZ4)
<a href="http://www.youtube.com/watch?feature=player_embedded&v=RAf4r4mEuZ4" target="_blank"><img src="http://img.youtube.com/vi/RAf4r4mEuZ4/0.jpg" 
 </a>
※画像をクリックするとYouTube動画につながります

コードの主要な部分は以下のとおり
def lorenz(var, t, p, r, b):は上記と同じ、そしてgen_sin(t)は以下のとおり

```py
def get_sin(t):
    p = 10
    b = 8/3
    r = 28
    #y0 = [1,0, 0]   #,0,0,0]
    y0 = [0,10,20]
    t1 = np.linspace(0, t-ts, 500)
    while len(t1)==0:
        t = time.time()
        t1 = np.linspace(0, t-ts, 500)
    try:
        sol = odeint(lorenz, y0, t1, args=(p, r, b))
    except:
        print('計算出来ませんでした！')
        print('t1= ',t1)
        print('y0= ',y0)
        sol = odeint(lorenz, y0, t1, args=(p, r, b))
    x, y, z = sol.T[0], sol.T[1], sol.T[2]
    return x,y,z
```
そして、これを受けるUnity側は[前回のもの](https://qiita.com/MuAuan/items/82dbdb673b799a615fd1)と以下の部分のみ変更した。

```py
float s = float.Parse(m_MyText.text.Substring(0,10));　
float s1 = float.Parse(m_MyText.text.Substring(startIndex: 11, length: 10));
float s2 = float.Parse(m_MyText.text.Substring(startIndex: 22, length: 10));
```
###まとめ
・Loranz方程式の数値解を求め、初期値の違いによる振舞をみた
・Lorenz方程式のカオスな振舞をOculusGoでVRした

・常微分方程式の解法とデータ授受が不安定であり、まとめて渡す方法を確立する必要がある