先ほどDuffing振動子のカオスな振舞をVRしたが、今回はカオス研究の発端となったLorenz方程式のカオスな振舞をVRしようと思う。
こちらは、Duffing振動子は一次元二階常微分方程式を解いたが、今回は一階の３次元常微分方程式を解く必要がある。
ちょっと目先が異なるのでodeintで解けるものやらと思っていたが、解けたので記載しておこうと思う。
やったこと
（１）Lorenz方程式の解法
（２）Lorenz方程式の振舞
（３）OculusGoでVRしてみる
（１）Lorenz方程式の解と振舞
コードは以下に置いた
・LSTM/duffing.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
あるいは、以下のようなものである。
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
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方程式の振舞
{\frac {dx}{dt}}=-px+py
{\frac {dy}{dt}}=-xz+rx-y
{\frac {dz}{dt}}=xy-bz
結果は、以下の参考の動画像の通りだが、。。やってみると結構面白い結果になった。ここでは以下代表的と思われる初期値等のセットを掲載する。上記のコードで比較的簡単に適当な初期値で計算できるので試せると思う。
【参考】
・ローレンツ方程式
・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_yz
plot_zx
そして時間変化は
plot_xt
plot_yt
plot_zt
一方、y0 = [0.1,0, 0.0]のとき
t = np.linspace(0, 20, 3001)とすると以下が得られる。
plot_xy
plot_yz
plot_zx
plot_xt
plot_yt
plot_zt
（３）OculusGoでVRしてみる
さて、いよいよVRで見えるかな？？
ほぼ前回と同じなので見えるはず。。
今回は、３パラメータの授受だけど基本同じ手法でやりました。
OculusGoでLorenz方程式のカオスな運動を見る♬
※画像をクリックするとYouTube動画につながります
コードの主要な部分は以下のとおり
def lorenz(var, t, p, r, b):は上記と同じ、そしてgen_sin(t)は以下のとおり
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側は前回のものと以下の部分のみ変更した。
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した
・常微分方程式の解法とデータ授受が不安定であり、まとめて渡す方法を確立する必要がある