【常微分方程式】Duffing振動子の解。。。♬

LSTMで株価フィッティングを調べていたら、先行記事があって、Duffing振動子のカオスの話。。。ということで、まずそこから入ることにしました。

【参考】
1.[【Python】QRNNでカオス時系列データ予測【Keras】]
(https://qiita.com/yukiB/items/681f68690ffabbf3e1e1)
2.Pythonでカオス・フラクタルを見よう！
3.Duffing振動子のカオス的運動のシミュレーション
4.scipy.integrate.odeint

やったこと

（１）常微分方程式の解
（２）Duffing振動子の解
（３）カオスということ
（４）LSTMでフィッティングするためのコード

（１）常微分方程式の解

これは、Runge-GittaとかRunge-Gitta－Gillとかで解けるけど、。。なんと参考4.のようにscipy.integrate.odeintとscipy.integrate.odeとかあるんですね。。。
ということで、参考2．や1．で使っているしかも使いやすそうなodeintというのを使ってみる。
参考4．の例によれば、通常の振り子の常微分方程式は以下のようなコードである。

"""
scipy.integrate.odeint
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html
The second order differential equation for the angle theta of a pendulum acted on by gravity with friction can be written:
theta''(t) + b*theta'(t) + c*sin(theta(t)) = 0
theta'(t) = omega(t)
omega'(t) = -b*omega(t) - c*sin(theta(t))
"""

from scipy.integrate import odeint, simps
import numpy as np
import matplotlib.pyplot as plt

def pend(y, t, b, c):
    theta, omega = y
    dydt = [omega, -b*omega - c*np.sin(theta)]
    return dydt

b = 0.025
c = 5.0
# 初期値y0=[theta0,omega0]
y0 = [np.pi - 0.1, 0.0]
# [0, 100]の範囲を1001の等間隔の時間tを生成
t = np.linspace(0, 100, 1001)
# odeint(関数,初期値,時間,関数の定数)
sol = odeint(pend, y0, t, args=(b, c))

# x=sol.T[0], p(dx/dt)=sol.T[1]
x, p = sol.T[0], sol.T[1]
# 以下グラフ化；theta(t)vs.tを書いたあと、omega(t)vs.tを上書き
plt.plot(x, p, ".", markersize=4)
plt.pause(3)
plt.close()
plt.plot(t, sol[:, 0], 'b', label='theta(t)')
plt.xlabel('t')
plt.grid()
plt.legend(loc='best')
plt.pause(3)
plt.plot(t, sol[:, 1], 'g', label='omega(t)')
plt.legend(loc='best')
plt.pause(3)
plt.close()

こんな簡単なコードで以下が得られます。

まず、x vs. tとp vs. tは以下のとおり、減衰振動の時間変化が得られます。

そして、p vs. xは以下のとおり減衰振動の位相空間での振る舞いが得られます。

（２）Duffing振動子の解

Duffing振動子の微分方程式は参考から以下のとおり、

m\frac{d^2x}{dt^2} = -\gamma \frac{dx}{dt} + 2ax - 4bx^3 + F_0\cos(\omega t + \delta)

ここで、各パラメータはそれぞれの参考で異なるけど、一応プログラム的には以下のとおり、書けます。
※これは参考2．のとおり

duffing.py

from scipy.integrate import odeint, simps
import numpy as np
import matplotlib.pyplot as plt

def duffing(var, t, gamma, a, b, F0, omega, delta):
    """
    var = [x, p]
    dx/dt = p
    dp/dt = -gamma*p + 2*a*x - 4*b*x**3 + F0*cos(omega*t + delta)
    """
    x_dot = var[1]
    p_dot = -gamma * var[1] + 2 * a * var[0] - 4 * b * var[0]**3 + F0 * np.cos(omega * t + delta)
    return np.array([x_dot, p_dot])

# parameter
F0, gamma, omega, delta = 10, 0.1, np.pi/3, 1.5*np.pi
a, b = 1/4, 1/2
# 参考3．のパラメーター
# F0, gamma, omega, delta = 2.0, 0.1, 2.4, 0
# a, b = 1/5, 1/2
# 参考2．の初期値
var, var_lin = [[0, 1]] * 2

# timescale　(0,20000)を2*np.pi/omega(=6)で区切って時間生成
t = np.arange(0, 20000, 2*np.pi/omega)
# t = np.linspace(0, 200, 3000)   #10000)　
# (0,100)を10000個の時間生成
t_lin = np.linspace(0, 100, 10000)

# solve　var, t の組に対して解を得る
var = odeint(duffing, var, t, args=(gamma, a, b, F0, omega, delta))
# solve var_lin, t_linの組に対して解を得る（描く絵の都合で二種類解を求めている）
# var_lin = odeint(duffing, var_lin, t_lin, args=(gamma, a, b, F0, omega, delta))
# 以下は参考3．の初期値x=0.5, dx/dt=0
var_lin = odeint(duffing, [0.5, 0], t_lin, args=(gamma, a, b, F0, omega, delta))

x, p = var.T[0], var.T[1]
x_lin, p_lin = var_lin.T[0], var_lin.T[1]

# 100個（つまりt=600）まで描画
plt.plot(t[:100], x[:100])
plt.pause(3)
# plt.close()
plt.plot(t[:100], p[:100])
plt.pause(3)
plt.savefig('plot_D_xpvst0.png', dpi=60)
plt.close()
plt.plot(x, p, ".", markersize=4)
plt.pause(3)
plt.savefig('plot_D_xp0.png', dpi=60)
plt.close()
# is chaotic?
plt.plot(t_lin, x_lin)
plt.pause(3)
# var_lin = odeint(duffing, [0.1, 1], t_lin, args=(gamma, a, b, F0, omega, delta))
# 以下は参考3．の初期値x=0.5001, dx/dt=0
var_lin = odeint(duffing, [0.5001, 0], t_lin, args=(gamma, a, b, F0, omega, delta))
x_lin, p_lin = var_lin.T[0], var_lin.T[1]
plt.plot(t_lin, x_lin)
plt.pause(3)
plt.savefig('plot_D_xxvst1.png', dpi=60)
plt.close()
plt.plot(x_lin, p_lin,".", markersize=4)
plt.pause(3)
plt.savefig('plot_D_xp1.png', dpi=60)
plt.close()

ということで、以下のようなグラフが得られます。
plot_D_xpvst0.png

plot_D_xp0.png

plot_D_xxvst1.png

plot_D_xp1.png

（３）カオスということ

「決定論的な連立常微分方程式が初期値鋭敏性を持つ」
カオスの定義あるいはカオスと呼ばれるものの特性とは、「非線形な決定論的力学系から発生する、初期値鋭敏性を持つ、有界な非周期軌道」
ということのようです。

ここで得られた結果plot_D_xp0.pngやplot_D_xp1.pngを見ると、上記の減衰振動の位相空間の振る舞いと比較すれば一目瞭然ですが、ちょっとした初期条件の違いが全く異なる地点に導かれてしまう。また、plot_D_xpvst0.pngなどを合わせて見れば確かに有界な非周期軌道であることも納得できます。

これが、株価と関係するかどうかはわかりません。
株価はもちろん決定論的な方程式に従ってはいません、しかしその変動も小さな揺らぎが拡大して異なる結果を導いているようにも見えるという意味ではカオス的なふるまいの時間変化をフィッティングして予測できれば、フィッティングできるかもと期待させてくれるということで価値があるように思います。
※それが出来たら気象にも応用できそうだし、逆に株価予測には気象のアンサンブル予報のような考え方も有効なのかもしれません

【参考】
・ジャパニーズ・アトラクタ
・ローレンツ方程式
・カオス理論

（４）LSTMでフィッティングするためのコード

ということで、回り道をしましたが、株価予測のためのコードの前に振り子の振る舞いをLSTMで予測するコードは、前回のgen_cosine_amp()関数の代わりに以下のようなものとなります。
※あとはこの常微分方程式をDuffing振動子のそれに変更するだけです

cos = np.zeros((5000, 1, 1))
def pend(y, t, b, c):
    theta, omega = y
    dydt = [omega, -b*omega - c*np.sin(theta)]
    return dydt
b = 0.025
c = 5.0
y0 = [np.pi - 0.1, 0.0]
t = np.linspace(0, 200, 5001)
sol = odeint(pend, y0, t, args=(b, c))
x, p = sol.T[0], sol.T[1]
plt.plot(x, p, ".", markersize=4)
plt.pause(3)
plt.savefig('plot_x-p_pendulum.png', dpi=60)
plt.close()

cos = np.zeros((5000, 1, 1))
for i in range(len(cos)):
    cos[i,0,0]=100*(x[i])/x[np.argmax(x)]
print('Generating Data...')
cos = cos   #gen_cosine_amp()

まとめ

・LSTMでDuffing振動子を予測するために、常微分方程式の解法を調査した
・振り子とDuffing振動子の振る舞いをグラフ化した
・振り子の解を利用してLSTMのInputとして利用する方法を確立した

・次回は、いよいよLSTMによる予測精度について記載したいと思う
・株価読み取りと予測を行う予定である