追記： この記事はTensorFlow のバージョン1系を使って書いたので，古い情報になっています． 最近の記事ではTensorFlow 2系を使ってJacobianの計算などをしていますが，この辺りは1と比べて随分簡単にできるようになっています．
#やりたいこと
Tensorflow を使って力学系の解析をします。今回はヤコビ行列を計算することにより，平衡点の安定性の解析、簡単な分岐の解析を行ってみます。
#環境
- Tensorflow 1.14
- scipy 1.2.1
#微分方程式系
前回1ローレンツ方程式の数値計算を行いました。
\frac{d \boldsymbol{x}(t)}{dt} = f(\boldsymbol{x}(t)) \tag{1}
とかける微分方程式を扱います。 $\boldsymbol{x}(t) \in \mathbb{R}^{d}$は$d$次元ベクトル、$f$は$\mathbb{R}^{d}$ から $\mathbb{R}^{d}$への写像（ベクトル場）です。
今回もローレンツ方程式を例として扱います。 $\boldsymbol{x}(t)=(x(t),y(t),z(t))^T$ として
\begin{equation}
\frac{d \boldsymbol{x}}{dt} =
\left( \begin{matrix}
f_x(\boldsymbol{x}) \\
f_y(\boldsymbol{x}) \\
f_z(\boldsymbol{x})
\end{matrix}\right) =
\left( \begin{matrix}
\sigma(y-x) \\
x(\rho -z) -y \\
xy - \beta z
\end{matrix}\right) \tag{2}
\end{equation}
$\sigma, \rho, \beta $ は定数のパラメータです。
この $f$ を計算するグラフを作る関数は以下です。
class Lorenz(object):
def __init__(self, sigma=10., beta=8 / 3., rho=28., **kwargs):
super().__init__(**kwargs)
self.sigma = sigma
self.beta = beta
self.rho = rho
def __call__(self, t, x):
dx_dt = self.sigma * (x[1] - x[0])
dy_dt = x[0] * (self.rho - x[2]) - x[1]
dz_dt = x[0] * x[1] - self.beta * x[2]
dX_dt = tf.stack([dx_dt, dy_dt, dz_dt])
return dX_dt
#ヤコビアン
(2)のようなベクトル値関数$f$に対する行列
D_{f} (\boldsymbol{x}) \equiv
\left( \begin{matrix}
\frac{\partial f_x}{\partial x}(\boldsymbol{x}) &
\frac{\partial f_x}{\partial y}(\boldsymbol{x}) &
\frac{\partial f_x}{\partial z}(\boldsymbol{x}) \\
\frac{\partial f_y}{\partial x}(\boldsymbol{x}) &
\frac{\partial f_y}{\partial y}(\boldsymbol{x}) &
\frac{\partial f_y}{\partial z}(\boldsymbol{x}) \\
\frac{\partial f_z}{\partial x}(\boldsymbol{x}) &
\frac{\partial f_z}{\partial y}(\boldsymbol{x}) &
\frac{\partial f_z}{\partial z}(\boldsymbol{x})
\end{matrix}\right)
を$f$のヤコビアン(ヤコビ行列、線形化行列とも）と呼びます。 一次元関数の微係数を一般化したものともいえます。元の点 $\boldsymbol{x}$ から微小な値 $\boldsymbol{\xi}$ ずれた点 $\boldsymbol{x}+\boldsymbol{\xi}$ についての時間発展は
\frac{d (\boldsymbol{x}+\boldsymbol{\xi})}{dt} = f(\boldsymbol{x}+\boldsymbol{\xi}) = f(\boldsymbol{x})+D_{f}(\boldsymbol{x}) \boldsymbol{\xi} + \mathcal{O}(\boldsymbol{|\xi}|^2)
とかけ、\boldsymbol{\xi} を無限小だとしてその時間発展を線形化すると線形化発展方程式
\frac{d \boldsymbol{\xi}}{dt} = D_{f}(\boldsymbol{x}(t)) \boldsymbol{\xi}
が得られます。
さて、$f$ からヤコビアンを求めるグラフを作成する関数は、TensorFlowでは次のように
することで実現できました。
def jacobian(t, f, x):
""" return jacobian matrix of f at x"""
n = x.shape[-1].value
fx = f(t, x)
if x.shape[-1].value != fx.shape[-1].value:
print('For calculating Jacobian matrix',
'dimensions of f(x) and x must be the same')
return
return tf.concat([tf.gradients(fx[i], x) for i in range(0, n)], 0)
ここで勾配を自動的に計算できるTensorFlowの恩恵を得ることができますね。この各要素の勾配を一度に計算できるやり方がなさそうだったので、次元の数だけ別個にgradientsを呼び出しています。
ヤコビアンの計算例
ローレンツ方程式の $f$ をつくり、それをもとにヤコビアン $D_f$ を作ります。
x = tf.placeholder(dtype=tf.float64, shape=(3), name='x')
f_lorenz = Lorenz()
Df = jacobian(None,f_lorenz,x)
placeholderのxに適当な値を入れて評価してみます。
sess = tf.Session()
sess.run(Df, feed_dict={x:[1,2,3]})
結果
array([[-10. , 10. , 0. ],
[ 25. , -1. , -1. ],
[ 2. , 1. , -2.66666667]])
手計算ではこの行列は、
\begin{equation}
D_{f} (\boldsymbol{x}) =
\left( \begin{matrix}
-\sigma &
\sigma &
0 \\
\rho-z &
-1 &
-x \\
y &
x &
-\beta
\end{matrix}\right) \tag{3}
\end{equation}
となるので、これに$\boldsymbol{x}=(1,2,3)$を入れると、ちゃんと計算できてることがわかると思います。
平衡点の安定性
$f(\boldsymbol{x}^{*})=\boldsymbol{0}$ となる点，つまり時間変化が0になる点 $\boldsymbol{x}^*$ を $f$の平衡点（不動点、固定点とも、）と呼びます。この点の上に解があると時間が進んでも解が動かない状態になります。
この点の近傍にある点が時間とともにこの点に近づいていくか、離れていくかを調べるのが安定性の解析です。（数学的な説明は厳密にはしてませんので、詳しくは教科書2などをみてください。）
ローレンツ方程式の自明な平衡点として原点 $(0,0,0)$ があるので、この点の安定性を調べてみましょう。固有値の計算は，TensorFlowのライブラリにある固有値を求める関数は自己随伴行列（実行列の場合対称行列）にしか使えないため、numpyの関数を使いました。
df0 = sess.run(Df, feed_dict={x:[0,0,0]}) #jacobian at origin
eigvals,eigvec = np.linalg.eig(df0)
print(eigvals)
結果
[-22.82772345 11.82772345 -2.66666667]
となり、固有値が３つ計算されました．実部が正の固有値が一つあるため，このときは原点は不安定であることがわかりました。 一般に，実部が正と負の固有値をともにもつ平衡点をサドルといいます。
分岐解析
上の結果はデフォルトのパラメータ $\sigma=10$, $\rho=28$, $\beta=8/3$ での結果です．このパラメータを変えて行った時に平衡点等の安定性が変わることがあります．平衡点の安定性が変わるなど，解の振る舞いが違う性質変わることを分岐といいます。
ここでは $\rho$ を変えていったときに原点の安定性がどう変わっていくかみてみましょう
以前はf_lorenzのrhoは浮動小数点型で与えていましたが，これをTensorFlowのplaceholder型としておきたいと思います．クラスのインスタンスを作るときに引数rhoにplaceholder型を渡すことで実現します．
x = tf.placeholder(dtype=tf.float64, shape=(3), name='x')
#bifurcation parameter
rho = tf.placeholder(dtype=tf.float64,
shape=(), name='rho')
f_lorenz = Lorenz(rho=rho)
Df = jacobian(None,f_lorenz,x)
この$\rho$ の値を0から2まで変化させていったときに固有値がどう変化するか見てみます．
sess=tf.Session()
sess.run(tf.initializers.global_variables())
rmin=0
rstep=0.02
rmax=2
li_rs=[]
li_eigvals=[]
for r in np.arange(rmin,rmax+rstep,rstep):
df = sess.run(Df, feed_dict={x:(0,0,0), rho:r})
eigvals,eigvecs = np.linalg.eig(df)
print("rho: {:.3g} lambda: {:.4g}".format(r,eigvals.real.max()))
li_rs.append(r)
li_eigvals.append(eigvals)
実行すると各$\rho$での固有値の実部のうち最大のものが表示されていきます．
途中を略しますが，
rho: 0.96 lambda: -0.03648
rho: 0.98 lambda: -0.01821
rho: 1 lambda: 0
rho: 1.02 lambda: 0.01815
rho: 1.04 lambda: 0.03624
となっており，$\rho=1$のとき最大の固有値が0となり$\rho>１$で原点が不安定になることがわかります．
グラフに書いてみます．
fig=plt.figure()
ax=fig.add_subplot(1,1,1)
ax.plot(li_rs,li_eigvals)
ax.grid()
ax.set_xlabel(r'$\rho$')
ax.set_ylabel(r'Re[$\lambda_i$]')
この結果は，解析的にも$ D_f(\boldsymbol{0})$ の行列式 $ \det D_f(\boldsymbol{0})$ を調べるなどして確かめられます．$\sigma>0, \beta>0$ の範囲では固有値が0となる条件は，$\rho=1$ であることがわかりますので，今回の数値計算の結果が正確であることが確認できます．
#まとめ
TensorFlowを使って微分方程式の $f$ を作っておくと，ヤコビアンを計算し，安定性の数値計算などを行うことが比較的短いコードで簡単にできることをみました．
#次回
今回は自明な平衡点があったので，その安定性をみましたが，平衡点がどこにあるかわからない場合，正確な位置が解析的にはわからない場合，どうすればよいでしょうか？どうにかして数値的に平衡点に近い近似解を見つけなければなりません．このために勾配が役に立つはずです．
次回はTensorFlowの最急降下法を使って力学系の平衡点を見つけるということをやってみたいと思います．
#code
githubに公開しています．
https://github.com/yymgch/ode_tf/blob/master/tf-dynamicalsystem-1-jacobian.ipynb
-
https://qiita.com/yymgt/items/4313d5ffe8d88486ed06
今回は引き続き ↩ -
Hirsch, Smale, Devaney 力学系入門 ―微分方程式からカオスまで―
平衡点 $\boldsymbol{x}^*$ の安定性を調べるにはこの点におけるヤコビアン $D_f(\boldsymbol{x}^*)$ の固有値を調べます．
そして全ての固有値の実部が負であるなら、平衡点は（漸近）安定であり，$x^*$ の近くの点は指数的な速さで $x^{*}$ に近づきます． ↩