📂動力学

ローレンツ・アトラクター

概要

ローレンツ方程式は、大気の対流を連立常微分方程式で表した数学的モデルである。

システム

$$ \begin{align*} {{dx} \over {dt}} =& - \sigma x + \sigma y \\ {{dy} \over {dt}} =& - xz + \rho x - y \\ {{dz} \over {dt}} =& xy - \beta z \end{align*} $$

変数

• $x(t)$: $t$ 時点の粒子の$x$座標を示す。
• $y(t)$: $t$ 時点の粒子の$y$座標を示す。
• $z(t)$: $t$ 時点の粒子の$z$座標を示す。

パラメータ

• $\sigma$: 粘性と熱伝導率に関するパラメータであるプラントル数だ。
• $\rho$: 流体の熱伝達方法に関するパラメータであり、レイリー数だ。
• $\beta$: ウィキペディア¹によると、次元に関連するパラメータらしい。教科書でもこれが何であるかを正確に説明しないことが多い。

---

初めのシミュレーションでは、これらのパラメータはローレンツが行ったように$\sigma = 10$、$\rho = 28$、$\displaystyle \beta = {{8} \over {3}}$ とすることが多い。レイリー数が $\rho \approx 24.74$ を超えると、ローレンツアトラクタが描く粒子の運動はカオスになる。非常に小さい初期条件の変化にも敏感に反応し、その予測は非常に困難だ。蝶が一回羽ばたくことで嵐が起こりうるということだ。

バイファーケーション

上の図はローレンツアトラクターのバイファーケーション図である。横軸は$25$から$325$まで少しずつ変わる$\rho$の値で、縦軸は十分に長い時間後の軌道で$z$方向の極大値をプロットしたものだ。$r= \rho$が小さいときは非常に不安定な動きを見せるが、一定レベル以上に大きくなると安定した軌道を持つようになることが分かる。次の説明を読むことをお勧めする：

上のグラフは、$\gamma=100$ の時の軌跡を描いたものである。通常、$\gamma$がある程度以上大きい場合、初めはカオス的な動きを見せるように思われるが、時間が経つにつれ周期的になることができる。

上のグラフは、最初の軌跡を取り除き、長い時間が経過した後の粒子の軌跡のみを示している。完全に周期的ではないが、ほぼ行動を予測できるレベルになった。

$xz$ のグラフを描くと、上のような形をしている。最初はカオス的に見えたが、時間が経つにつれて軌道が安定していき、極大値は３つほどにはっきりと現れる。実際にバイファーケーション図を見ても同じ解釈を出すことができる。

$\gamma = 28$ のときは、十分に長い時間が経っても次のようにカオス的な動きを見せる。バイファーケーション図では、多くの点が散らばって厚くなった形状を成す。

一方、$\gamma = 300$のときは、次のように明確に現れる２つの極大値を持つ。バイファーケーション図では、一つの$\gamma$値に対して２点が対応することで、２つの線が現れるようになる。

コード

MATLAB

以下はMATLABを使ってローレンツアトラクターの軌跡を描くコードである。

sigma = 10;
beta = 8/3;
rho = 28;
f = @(t,v) [-sigma*v(1) + sigma*v(2); rho*v(1) - v(2) - v(1)*v(3); -beta*v(3) + v(1)*v(2)];
[t,v] = ode45(f,[0 10000],[1 1 1]);     % Runge-Kutta 4th/5th order ODE solver
x=v(:,1); x=x(50000:end);
y=v(:,2); y=y(50000:end);
z=v(:,3); z=z(50000:end);
plot3(x,y,z)
%plot(x,z)

Fortran

以下はFortranを使ってバイファーケーション図のソースを作り、OUTPUT.csvファイルとして返すコードだ。ODEを解くために使用された方法はRK4だ。

program LorenzBifurcation
  implicit none
  double precision :: sigma=10.d0, rho=28.d0, beta=8.d0/3.d0
  Integer, parameter :: D = 3
  double precision :: v(D)=1.d0, h=0.001d0, z1, z2, z3
  integer :: i, j, k, preitr=10**4-1000, endstep=10**4
  character :: answer
  abstract interface
     function func (v)
       Integer, parameter :: D = 3
       double precision :: v(D), func(D)
     end function func
  end interface
  open(UNIT=34,FILE="OUTPUT.csv",ACTION="WRITE")
  rho=5.d0
  h=0.01d0
  do while(rho<325)
    print '(A,f6.4,A)', "Now", (rho-5.d0)/320.d0, "%"
    do i=1, preitr
       v = RK4(v,f)
    end do
    z1 = preitr
    z2 = preitr
    z3 = preitr
    do i=preitr+1, endstep
       v = RK4(v,f)
       z1 = z2
       z2 = z3
       z3 = v(3)
       if(z1<z2 .and. z2>z3) then
         write (34,*) rho, ",", z2
       end if
    end do
    v=1.d0
    rho=rho+0.1d0
  end do
contains
  function RK4(v,f)
    procedure (func) :: f
    double precision :: v(D), V1(D), V2(D), V3(D), V4(D), RK4(D)
    V1 = f(v)
    V2 = f(v + (h/2)*V1)
    V3 = f(v + (h/2)*V2)
    V4 = f(v + h*V3)
    RK4 = v + (h/6)*(V1 + 2*V2 + 2*V3 + V4)
  end function RK4
  
  function f(v)
    double precision :: v(D), f(D)
    f(1) = sigma*(v(2)-v(1))
    f(2) = v(1)*(rho-v(3))-v(2)
    f(3) = v(1)*v(2) - beta*v(3)
  end function f
end program LorenzBifurcation

以下は、上で得られたOUTPUT.csvファイルを使ってバイファーケーション図を描くMATLABのコードだ。

[file,path] = uigetfile('*.csv');
data = csvread([path,file]);
 
sz = 0.6;
scatter(data(:,1),data(:,2),sz,'.');

Julia

以下はJuliaでローレンツアトラクターを見つけて描くコードである。

using DifferentialEquations
using Plots

function parameterized_lorenz!(du,u,p,t)
    x,y,z = u
    σ,ρ,β = p
    du[1] = dx = σ*(y-x)
    du[2] = dy = x*(ρ-z) - y
    du[3] = dz = x*y - β*z
end
u0 = [1.0,0.0,0.0]
tspan = (0.0,100.0)
p = [10.0,28.0,8/3]
prob = ODEProblem(parameterized_lorenz!,u0,tspan,p)
solution = solve(prob)
plot(solution, vars = (1,2,3))