Pythonでローレンツアトラクタをかいたよ
Pythonの練習のためにローレンツアトラクタを描きました。
ローレンツアトラクタ(Lorenz attractor)とは? >> ローレンツ方程式 - Wikipedia
カオスの教科書の一番最初に登場するやつです。めっちゃ単純な方程式なのにパラメータによってめちゃくちゃ解の挙動が変わるところが面白いです。
前にC++で描いたことがあるんですけど、PythonだとMatplotlibとかいう描画ライブラリを使えば数値計算とグラフ描画を同時にやってくれるのがとてもありがたいですね。
パラメータ、初期値でRungeKutta法を使って解いたときの結果です。
Pythonのコードです。まだCっぽい書き方になっちゃってます。
import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # const p = 10 r = 28 b = 8/3 dt = 0.01 t_0 = 0 t_1 = 50 X_0 = np.array([1, 1, 1]) # function def RungeKutta(t, X): k_1 = LorenzEquation(t, X) k_2 = LorenzEquation(t + dt/2, X + k_1*dt/2) k_3 = LorenzEquation(t + dt/2, X + k_2*dt/2) k_4 = LorenzEquation(t + dt, X + k_3*dt) X_next = X + dt/6*(k_1 + 2*k_2 + 2*k_3 + k_4) return X_next def LorenzEquation(t, X): x = X[0] y = X[1] z = X[2] return np.array([-p*x + p*y, -x*z + r*x - y, x*y - b*z]) # main process t = t_0 X = X_0 data = np.r_[X] while t < t_1: X = RungeKutta(t, X) t += dt data = np.c_[data, X] print(data) fig = plt.figure() ax = Axes3D(fig) ax.plot(data[0,:], data[1,:], data[2,:]) plt.show()
比較参考までにEigenライブラリを使ったC++のときのコードです。(Lorenzのつづりが間違ってる・・・)
#include <iostream> #include <fstream> #include <Eigen/Core> using namespace Eigen; Vector3f RungeKutta(float t, Vector3f X); Vector3f LorentzEquation(float t, Vector3f X); const float p = 10; const float r = 28; const float b = 8/3; const float dt = 0.01; const float t_0 = 0; const float t_1 = 10; const Vector3f X_0(1, 1, 1); int main() { float t = t_0; Vector3f X = X_0; std::ofstream ofs; ofs.open("lorentz_attractor.txt", std::ios::out); ofs << t << " " << X.transpose() << std::endl; do { X = RungeKutta(t, X); t += dt; ofs << t << " " << X.transpose() << std::endl; } while(t < t_1); return 0; } Vector3f RungeKutta(float t, Vector3f X) { Vector3f k_1 = LorentzEquation(t, X); Vector3f k_2 = LorentzEquation(t + dt/2, X + dt/2*k_1); Vector3f k_3 = LorentzEquation(t + dt/2, X + dt/2*k_2); Vector3f k_4 = LorentzEquation(t + dt, X + dt*k_3); Vector3f X_next = X + dt/6*(k_1 + 2*k_2 + 2*k_3 + k_4); return X_next; } Vector3f LorentzEquation(float t, Vector3f X) { float x = X[0]; float y = X[1]; float z = X[2]; Vector3f val(-p*x + p*y, -x*z + r*x - y, x*y - b*z); return val; }
今度は時刻によってパラメータを動的に変えたときの描画とか他のアトラクタも描いてみます。
attractorって日本語でなんて訳すのかな。