ローレンツ方程式：
ローレンツ方程式は、
dx/dt = -p*x + p*y
dy/dt = -x*z + r*x - y
dz/dt = x*y - b*z
であり、これら三つの微分方程式に任意の初期値x,y,zを与えると（x=0, y=0, z=0以外）、係数p=10, r=28, b=8/3の時にカオス的な振る舞いをするというもの。通常の数式なら発散/収束/周期的になりますが、そうはならないというところが面白い。また初期値がほんの少し違うだけで結果が大きく変わるという初期値鋭敏性も持つ。
初期値x=1, y=1, z=1を与えて、1ステップ0.001で50000ステップ分計算させたもの。
コードは以下。
import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D
dt = 0.001 T = 50 p, r, b = 10, 28, 8/3 print('steps:', int(T/dt)) def lorenz_euler(V): x, y, z = V[0], V[1], V[2] X, Y, Z = [], [], [] t = 0 while t <= T: dx = (-p*x + p*y) * dt dy = (-x*z + r*x - y) * dt dz = (x*y - b*z) * dt x += dx y += dy z += dz X.append(x) Y.append(y) Z.append(z) t += dt return X, Y, Z V0 = [1.0, 1.0, 1.0] X, Y, Z = lorenz_euler(V0) fig = plt.figure(figsize=(8, 8)) ax = fig.add_subplot(111, projection='3d') g = 25 ax.set_xlim(-g,g) ax.set_ylim(-g,g) ax.set_zlim(0, g*2) ax.plot(X, Y, Z, color='blue', alpha=0.6) plt.show()
とりあえず、xの式に注目すると、
dx/dt = -p*x + p*y
の式を変形して、
dx = (-p*x + p*y) * dt
にすれば、dxはxにおける1ステップ分の変化量なので、whileループ内で変化量dxをxに加算していき、設定した限度を越えればwhileループを打ち切り。
x += dx
という感じで、yやzについても同様に処理し、プログラミング的には簡単な手続き。
しかし調べてみると、これは「オイラー法」という計算方法であり、あまり精度が保てないらしい。精度を出すには1ステップの刻み幅をより細かくすれば単純に向上しますが、そうすると全体のステップ数も増え計算量が多くなります。そのための工夫としてルンゲクッタ法というのがあります。以下では、各計算方法の比較をしてみます。
オイラー法：
まず、簡単な例としてオイラー法から。
f(x) = x*sin(x)
という式で試してみます。
f(x)を微分すると、
f'(x) = sin(x) + x*cos(x)
になります。
xの変化量dxは、
dx = f'(x)*dt = (sin(x) + x*cos(x))*dt
となり、
x += dx
で、xの値を更新してループ処理させます。
以下コード。
h = 0.1 T = np.pi * 2 print('steps:', int(T/h)) def f(x): return x * np.sin(x) X0 = np.arange(0, T, h) Y0 = f(X0) def df(x): # derivative of x*sin(x) : (x*sin(x))' = sin(x) + x*cos(x) return np.sin(x) + x * np.cos(x) x = 0 y = 0 X = [x] Y = [y] while x <= T - h: dx = h dy = df(x) * h x += dx y += dy X.append(x) Y.append(y) fig = plt.figure(figsize=(8, 4)) ax = fig.add_subplot(111) ax.grid() ax.plot(X0, Y0, label='True ', alpha=0.8, color='C0') ax.plot(X, Y, label='Euler', alpha=0.8, color='C1') ax.legend()
1ステップの刻み幅は0.1。合計で62ステップ。X, Yが積分したオイラー法の結果で、X0, Y0は積分を使わないで求めた結果。
わずかですが誤差があります。刻み幅をもっと細かくすれば誤差も減りますが、ローレンツアトラクタの場合は初期値鋭敏性もあるためこの誤差は致命的かもしれません。
y値最大誤差：0.3043499782
２次のルンゲクッタ法：
この方法は改良オイラー法とも呼ばれているようで、その名の通りオイラー法よりは誤差が少ない。オイラー法では刻み幅とその時のy値（高さ）を掛け合わせているため、細長い長方形を積算していますが、２次のルンゲクッタ法では長方形上端を斜め（変化量の傾き）にした台形で計算します。その分精度向上するというもの。
h = 0.1 T = np.pi * 2 print('steps:', int(T/h)) def f(x): return x * np.sin(x) X0 = np.arange(0, T, h) Y0 = f(X0) def df(x): # derivative of sin(x) return np.sin(x) + x * np.cos(x) x = 0 y = 0 X2 = [x] Y2 = [y] while x <= T - h: k1 = df(x) k2 = df(x + h) x += h y += h / 2 * (k1 + k2) X2.append(x) Y2.append(y) fig = plt.figure(figsize=(8, 4)) ax = fig.add_subplot(111) ax.grid() ax.plot(X0, Y0, label='True ', alpha=0.8) ax.plot(X2, Y2, label='RK2', alpha=0.8) ax.legend()
k1とk2が台形の異なる二辺の長さ。あとは台形面積を求める計算。
今回はdtの代わりにhを刻み幅にしていますが同じ0.1です。誤差は0(h^3)程度。
かなり誤差はなくなりました。
y値最大誤差：0.0042499318
4次のルンゲクッタ法：
さらに精度をあげた方法。通常この4次のルンゲクッタ法くらいの精度が必要らしい。先ほどの2次ではk1、k2の２点でしたが、今回はk1〜k4の4点をとり、その平均を求めます。
h = 0.1 T = np.pi * 2 print('steps:', int(T/h)) def f(x): return x * np.sin(x) X0 = np.arange(0, T, h) Y0 = f(X0) def df(x): # derivative of x*sin(x) : (x*sin(x))' = sin(x) + x*cos(x) return np.sin(x) + x * np.cos(x) x = 0 y = 0 X4 = [x] Y4 = [y] while x < T - h: k1 = df(x) * h k2 = df(x + h/2) * h k3 = df(x + h/2) * h k4 = df(x + h) * h x += h y += (k1 + 2*k2 + 2*k3 + k4) / 6 X4.append(x) Y4.append(y) fig = plt.figure(figsize=(8, 4)) ax = fig.add_subplot(111) ax.grid() ax.plot(X0, Y0, label='True ', alpha=0.8) ax.plot(X4, Y4, label='RK4', alpha=0.8) ax.legend()テイラー展開により導出できるようですが、そのまま公式を覚えた方がよさそう。式中では刻み幅の両端点と中点をとって、1:2:2:1の割合で平均近似しているようです。そうすると、曲線上の中点の値df(x+h/2)と刻み幅hの両端を直線で結んだ時の中間平均値：(df(x) + df(x+h))/2には差分があるので、その２点の2/3の地点が計算で求まるようです。
グラフを見ただけでは違いはわかりにくいですが、数値的に誤差を求めてみるとさらに精度が上がっているようです。理論上、誤差は0(h^5)。
y値最大誤差：0.0000003046
自前でルンゲクッタ法をコーディングするかわりに、scipyのライブラリを使ってもいいかもしれません。
scipy.integrate.solve_ivp （methodでRK45など選択可）
scipy.integrate.odeint （古いタイプ？）
scipy.integrate.RK4 （使いにくい）
solve_ivpでは計算法としてRK23やRK45も選択可能。t_evalで刻み幅を配列で指定。さらにオプションのrtol, atolで精度を指定することもでできるようなので、アルゴリズム的な工夫はほぼ必要なし。デフォルトではrtol=1e-3、atol=1e-6になっており、そのままでも大丈夫そう。
RK2、RK4、Odeint、Solve_ivpの比較：
ルンゲクッタ法（RK2、RK4）、ならびにライブラリのOdeint、Solver_ivp（RK45）をローレンツアトラクで比較してみることに（とは言っても、初期値鋭敏性のあるローレンツアトラクタなので、精度の比較はしにくいですが）。
さらに、3D/2D表示、XYZの各値の変化を1D表示するファンクションも追加しておきました。
import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy.integrate import odeint from scipy.integrate import solve_ivp %matplotlib inline #%matplotlib notebook def df(t, xyz): x, y, z = xyz return np.array([-p*x + p*y, -x*z + r*x - y, x*y - b*z]) def lorenz_RK2(V): xyz = V t = 0 X, Y, Z = [], [], [] while t <= T: k1 = df(t, xyz) k2 = df(t+h, xyz + h*k1) xyz += (k1 + k2) * h / 2 X.append(xyz[0]) Y.append(xyz[1]) Z.append(xyz[2]) t += h return X, Y, Z def lorenz_RK4(V): xyz = V t = 0 X, Y, Z = [], [], [] while t <= T: k1 = df(t, xyz) k2 = df(t+h/2, xyz + h*k1/2) k3 = df(t+h/2, xyz + h*k2/2) k4 = df(t+h, xyz + h*k3) xyz += (k1 + 2*k2 + 2*k3 + k4) * h / 6 X.append(xyz[0]) Y.append(xyz[1]) Z.append(xyz[2]) t += h return X, Y, Z def lorenz_Odeint(V): xyz = V t = np.arange(0, T, h) X, Y, Z = odeint(df, xyz, t, tfirst=True).T #X, Y, Z = odeint(f2, xyz, t, rtol=1e-9, atol=1e-9).T return X, Y, Z def lorenz_ivp(V): xyz_0 = np.array(V) t = np.arange(0, T, h) res = solve_ivp(df, (0, T), xyz_0, method='RK45', t_eval=t)#, rtol=1e-7, atol=1e-9) X = res.y[0] Y = res.y[1] Z = res.y[2] print(X[-1], Y[-1], Z[-1]) return X, Y, Z def plot3D(methods, labels): from matplotlib.colors import BASE_COLORS colorlist = list(BASE_COLORS.keys()) fig = plt.figure(figsize=(10,10)) ax = fig.add_subplot(111, projection='3d') g = 30 ax.set_xlim(-g, g) ax.set_ylim(-g, g) ax.set_zlim(0, 2*g) ax.xaxis.pane.set_facecolor('white') ax.yaxis.pane.set_facecolor('white') ax.zaxis.pane.set_facecolor('white') for i, XYZ in enumerate(methods): X = XYZ[0] Y = XYZ[1] Z = XYZ[2] #color = 'C' + str(i) color = colorlist[i] ax.plot(X, Y, Z, color=color, alpha=0.5, label=labels[i]) # plot method ax.scatter(X[0], Y[0], Z[0], color=color, s=60, marker='x') # start point ax.text(X[0], Y[0], Z[0], s='START ', color=color, verticalalignment='top', horizontalalignment='right') ax.scatter(X[-1], Y[-1], Z[-1], color=color, s=30, marker='o') # end point plt.legend() plt.show() def plot2D(methods, labels, axis): from matplotlib.colors import BASE_COLORS colorlist = list(BASE_COLORS.keys()) fig = plt.figure(figsize=(8,8)) ax = fig.add_subplot(111) ax.axis([-30, 30, 0, 60], 'equal') if axis == 'XY': ax.axis([-30, 30, -30, 30], 'equal') ax.grid() for i, XYZ in enumerate(methods): if axis == 'XY': A = XYZ[0] B = XYZ[1] elif axis == 'XZ': A = XYZ[0] B = XYZ[2] else: A = XYZ[1] B = XYZ[2] color = colorlist[i] ax.plot(A, B, color=color, alpha=0.5, label=labels[i]) # plot line ax.scatter(A[0], B[0], color=color, s=60, marker='x') # start point ax.text(A[0], B[0], s='START ', color=color, verticalalignment='center', horizontalalignment='right') ax.scatter(A[-1], B[-1], color=color, s=30, marker='o') # end point plt.legend() plt.show() print(axis + ' coordinates') def plot1D(methods, labels, axis): plt.figure(figsize=(14, 3)) for i, XYZ in enumerate(methods): if axis == 'X': A = XYZ[0] elif axis == 'Y': A = XYZ[1] else: A = XYZ[2] plt.plot(list(range(len(A))), A, color='C'+str(i), alpha=0.6, label=labels[i]) plt.legend() plt.show()
実行用のコードは以下。
h = 1e-3 # dt T = 30 # total time p, r, b = 10, 28, 8/3 print('steps:', int(T/h)) xyz_0 = [10.0, 10.0, 10.0] # initial values for the methods methods = [lorenz_RK2(xyz_0), lorenz_RK4(xyz_0), lorenz_Odeint(xyz_0), lorenz_ivp(xyz_0)] labels = ['RK2', 'RK4', 'ODE', 'IVP'] plot3D(methods, labels) # plot2D(methods, labels, 'XY') # plot2D(methods, labels, 'XZ') plot2D(methods, labels, 'YZ') plot1D(methods, labels, 'X')
表示結果。
YZ座標。この方向からみるとほぼ線対称の形になります。
6000ステップくらいまでは一緒ですが、その後は次第にずれてきています。単体だけ見ると、規則的に見える部分もありますが、やはり非周期的になっています。
初期値鋭敏性：
計算方法は同じで、初期値をほんの少しだけ変えて実験してみました。
初期値A：X=10.0、Y=10.0、Z=10.0
初期値B：X=10.001、Y=10.0、Z=10.0
BのXの初期値だけ0.001違います。あとは同じ条件で50000ステップ。
以下が実行コード。
h = 1e-3 T = 50 p, r, b = 10, 28, 8/3 print('steps:', int(T/h)) xyz_0 = [10.0, 10.0, 10.0] xyz_1 = [10.001, 10.0, 10.0] methods = [lorenz_RK4(xyz_0), lorenz_RK4(xyz_1)] labels = ['RK4:'+str(xyz_0), 'RK4:'+str(xyz_1)] plot3D(methods, labels) # plot2D(methods, labels, 'XY') # plot2D(methods, labels, 'XZ') # plot2D(methods, labels, 'YZ') plot1D(methods, labels, 'X') plot1D(methods, labels, 'Z')
結果表示。
10000ステップを超えたあたりから差がでてきます。ずれるというより、異なったパターンになってきています。
完璧なランダムとは言えませんが、数多くのステップ数を繰り返しても周期的にならないという部分が興味深い。
分布について：
ランダムの指標として、値全体の分布もとってみました。
基本的には中央に寄っていますが、正規分布とも違う。
まとめ：
元々はカオス的な振る舞いをするローレンツアトラクタを利用して乱数生成できないかということを考えておりましたが、計算方法であるルンゲクッタ法やライブラリについての実験となりました。これはこれでいろいろ勉強になります。精度を高める計算方法にはいくつかあるようですが、PID制御もそのうちの一つという感じです。ということで、引き続き調査してみたいと思います。
