語り得ぬことについては、沈黙するほかない

世界は事実の総体であり、ものの総体ではない

Pythonでローレンツアトラクタをかいたよ

Pythonの練習のためにローレンツアトラクタを描きました。

ローレンツアトラクタ(Lorenz attractor)とは? >> ローレンツ方程式 - Wikipedia

カオスの教科書の一番最初に登場するやつです。めっちゃ単純な方程式なのにパラメータによってめちゃくちゃ解の挙動が変わるところが面白いです。

前にC++で描いたことがあるんですけど、PythonだとMatplotlibとかいう描画ライブラリを使えば数値計算とグラフ描画を同時にやってくれるのがとてもありがたいですね。

パラメータ \displaystyle (p, r, b) = (10, 28, 8/3)、初期値 \displaystyle (x_0, y_0, z_0) = (1, 1, 1)でRungeKutta法を使って解いたときの結果です。

f:id:hoshimure-47:20170801072843p:plain


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;  
}


今度は時刻 tによってパラメータを動的に変えたときの描画とか他のアトラクタも描いてみます。

attractorって日本語でなんて訳すのかな。