Lorenz方程式

記事の内容


  • Lorenz方程式
    • 固有値による解析
    • アトラクタの観察
    • 時系列のプロット
    • Lorenz写像
    • カオス的な振る舞い
  • ソースコード
    • 共通のコード
    • その他のプロット

今回はLorenz方程式を定性的に解析します.

Lorenz方程式

Lorenz方程式は次式の系です.

\begin{align} \dot{x} &= \sigma (y-x)\\ \dot{y} &= rx-y-xz\\ \dot{z} &= xy-bz \end{align}

以下, 3つのパラメータの値を次のように固定します.

\[ \sigma = 10,\quad r=28,\quad b=\frac{8}{3} \]

固有値による解析

この方程式の不動点は3つです.

\[ (0,0,0),\quad (6\sqrt{2}, 6\sqrt{2}, 27),\quad (-6\sqrt{2}, -6\sqrt{2}, 27) \] それぞれの不動点の周りで系を線形化し, 固有値を数値的に求めると次のようになります. \[ \lambda_1 = -22.828,\quad \lambda_2 = -0.375,\quad \lambda_3 = 11.828\\ \] \[ \lambda_1 = -11.32 - 5.687i,\quad \lambda_2 = -11.32 + 5.687i,\quad \lambda_3 = 8.973\ \] \[ \lambda_1 = -13.174, \quad \lambda_2 = -0.246 - 8.178i, \quad \lambda_3 = -0.246 + 8.178i\\ \]

よって, 3つ目の不動点を除いて, 不安定な不動点です.

アトラクタの観察

Lorenz方程式の軌道を描いてみます. 軌道の計算には4段4次Runge-Kuttaを用いました. 反復回数は10,000回, 刻み幅は0.01, 初期値は(0.0, 1.0, 0.0)としました. 以下は計算結果をプロットしたものです. 有名なストレンジアトラクタが観察できます. 軌道は周期的ではなく, カオス的な振る舞いを示します.

【コード3の実行結果】

時系列のプロット

時系列変化を調べます. 以下に, xとyの変化を示します. zに関しては次節で示します.

【コード4の実行結果】

xとyの時系列変化は一致しているように見えます. 実際には少しずれています. また, 正の部分と負の部分を行き来しています.

Lorenz写像

zの極大点の変化を追ってみます. コード6,7の結果の図は, 極大点(オレンジ色)の位置を表しています. コード5の結果は, zの極大点のみを取り出し, 次の極大点との組を点列としてプロットしたものです. この写像 \[ T(z_n) = z_{n+1} \] をLorenz写像といいます. 図を見ると, 直線的に右肩上がりになっている部分があります. この部分は恒等写像に見えるので, Lorenz写像の不動点が存在を示唆します.

【コード5の実行結果】
【コード6の実行結果】
【コード7の実行結果】

ここまでパラメータrを固定していました. rを動かしたときの軌道を観察すると, カオス的な振る舞いが観察できます. 次の図は, 各rの値(横軸)に対して, zの極大値(縦軸)をプロットしたものです. 白く抜けている部分があります. 縦に黒点が散らばっているほど, 軌道は極大点周りで様々な値をとります. 逆に, 白く抜けている部分では軌道がほぼ同じ点を通ります. これはLorenz写像がいくつかの不動点や周期点を持つことに対応します. 以下の図でこれらの事実を観察します. 下図から, r=150付近とr=200以降は軌道が重なっていると予想されます. なお, r=340付近のように, 1点に潰れている部分と, 2つに分かれている場合があります. 2つに分かれている場合には, 周期的に2つの値をとり続けます.

【コード8の実行結果】

以下は, r=340に固定した場合の種々のプロットです. 下のLorenz写像は不動点を持ちます.

【コード9の実行結果】

zの軌道も同じ値をとり続けます.

【コード9の実行結果】

軌道がほぼ重なっています.

【コード9の実行結果】

カオス的な振る舞い

Lorenz方程式はカオス的な振る舞いを示します. これは, 初期値鋭敏性からも分かります. 初期値を少しだけずらした2つの軌道を追ってみましょう. 以下の図は, 2つの初期値 \[ (0.0, 1.0, 0.0),\quad (10^{-10}, 1.0+10^{-10}, 10^{-10}) \] から出発した2つの軌道の差のノルムをプロットしたものです. 時間発展に伴って, 差はどんどん大きくなっていきます. 対数プロットで直線的に変化しているので, 指数関数的に差が広がります.

【コード10の実行結果】

ソースコード

共通のコード

【Juliaコード1; インポート】
using Plots
using LinearAlgebra
【Juliaコード2; 関数などの用意】
#パラメータ
σ = 10.0;
r = 28.0;
b = 8/3;

#RK法のパラメータ
n = 10000; #反復回数
h = 0.01; #刻み幅
x₀ = [0.0; 1.0; 0.0]; #初期値

#右辺関数
f(x, σ, r, b) = [σ*(x[2]-x[1]); r*x[1]-x[2]-x[1]*x[3]; x[1]*x[2]-b*x[3]];

#4段4次RK法の1ステップ
function one_step_RK(x_old::Array{Float64,1}, h::Float64, f, σ::Float64, r::Float64, b::Float64)
    #kの更新
    k₁ = f(x_old, σ, r, b);
    k₂ = f(x_old.+h*k₁/2, σ, r, b);
    k₃ = f(x_old.+h*k₂/2, σ, r, b);
    k₄ = f(x_old.+h*k₃, σ, r, b);

    #xの更新
    x_new = x_old .+ h*(k₁.+2*k₂.+2*k₃.+k₄)/6;
    return x_new
end

#RK法によるシミュレーション
function myRK(n::Int64, x₀::Array{Float64,1}, f, σ::Float64, r::Float64, b::Float64)
    #変数の初期化
    xvec = zeros(length(x₀),n);
    xvec[:,1] = x₀;

    #反復計算
    for k in 1:n-1
        x_old = xvec[:,k];

        #qとvの更新
        x_new = one_step_RK(x_old, h, f, σ, r, b);

        #更新値の保存
        xvec[:,k+1] = x_new;
    end
    return xvec
end

#極大点のインデックスを探す
function find_maximal_index(array)
    idx = [];
    for i in 2:length(array)-1
        #極大なら追加
        if (array[i-1] < array[i]) &(array[i+1] < array[i])
            append!(idx, i)
        end
    end
    return idx
end

その他のプロット

【Juliaコード3; アトラクタ】
#RK法
xvec = myRK(n, x₀, f, σ, r, b);
Lorenz_plot = plot(xvec[1,:], xvec[2,:], xvec[3,:], labels=false, title="Lorenz equation", xlabel="x", ylabel="y", zlabel="z");
#savefig(Lorenz_plot, "lorenz-figs/fig1.png");
【Juliaコード4; xとyの時系列】
#時系列のプロット
time_series_xy = plot(xvec[1,:], labels="x", size=(800,400), xlabel="time", ylabel="x or y");
plot!(xvec[2,:], labels="y");
#savefig(time_series_xy, "lorenz-figs/fig2.png");
【Juliaコード5; Lorenz写像】
#次項との関係
idx = find_maximal_index(xvec[3,:]);
maximal_z_array = xvec[3,idx];
lorenz_map_plot = scatter(maximal_z_array[1:end-1], maximal_z_array[2:end], labels="maximal point", xlabel="n", ylabel="n+1");
#savefig(lorenz_map_plot, "lorenz-figs/fig3.png");
【Juliaコード6; zの時系列】
#時系列のプロット
time_series_z = plot(xvec[3,:], labels=false, size=(800,400), xlabel="time", ylabel="z");
scatter!(idx, xvec[3,idx], markerstrokewidth=0, labels=false);
savefig(time_series_z, "lorenz-figs/fig4.png");
【Juliaコード7; Poincaré断面とLorenzアトラクタ】
#RK法
xvec = myRK(n, x₀, f, σ, r, b);
x_maximal = xvec[1,idx];
y_maximal = xvec[2,idx];
z_maximal = xvec[3,idx];
Lorenz_plot = plot(xvec[1,:], xvec[2,:], xvec[3,:], labels=false, title="Lorenz equation", xlabel="x", ylabel="y", zlabel="z");
scatter!(x_maximal, y_maximal, z_maximal, labels=false, markerstrokewidth=0);
savefig(Lorenz_plot, "lorenz-figs/fig5.png");
【Juliaコード8; パラメータを変化させる】
#rの最大値
rmax = 350.;

#プロット
maximal_plot = plot(legend=false);

#各rにおける極大値の計算
for r₀ in 1:0.5:rmax
    #軌道を計算する
    xvec = myRK(n, x₀, f, σ, r₀, b);

    #極大となるインデックス
    idx = find_maximal_index(xvec[3,:]);

    #保存
    scatter!(r₀*ones(length(idx)), xvec[3,idx], markersize=1, color=:black);
end

#保存
savefig(maximal_plot, "lorenz-figs/fig6.png");
【Juliaコード9; 閉軌道のプロット】
#rの値を固定する
r₀ = 160.0;

#RK法により軌道を計算
xvec = myRK(n, x₀, f, σ, r₀, b);

#Lorenz写像のプロット
idx = find_maximal_index(xvec[3,:]);
maximal_z_array = xvec[3,idx];
lorenz_map_plot = scatter(maximal_z_array[1:end-1], maximal_z_array[2:end], labels="maximal point", xlabel="zn", ylabel="zn+1");

#時系列のプロット, 極大点をオレンジで示す
time_series_z = plot(xvec[3,:], labels=false, size=(800,400), xlabel="time", ylabel="z");
scatter!(idx, xvec[3,idx], markerstrokewidth=0, labels=false);
x_maximal = xvec[1,idx];
y_maximal = xvec[2,idx];
z_maximal = xvec[3,idx];
Lorenz_plot = plot(xvec[1,:], xvec[2,:], xvec[3,:], labels=false, title="Lorenz equation", xlabel="x", ylabel="y", zlabel="z");
scatter!(x_maximal, y_maximal, z_maximal, labels=false, markerstrokewidth=0);

#保存
savefig(lorenz_map_plot, "lorenz-figs/fig7.png");
savefig(time_series_z, "lorenz-figs/fig8.png");
savefig(Lorenz_plot, "lorenz-figs/fig9.png");
【Juliaコード10; カオス的な振る舞い】
#初期値を少しずらす
x₀_new = x₀ .+ 1e-10*ones(3);

#それぞれの初期値で計算する
xvec = myRK(n, x₀, f, σ, r, b);
xvec_new = myRK(n, x₀_new, f, σ, r, b);

#軌道の差のノルムを時系列プロット
orb_diff = [norm(xvec[i].-xvec_new[i]) for i in 1:n];
diff_plot = plot(orb_diff, linewidth=0.1, yaxis=:log, xlabel="time", ylabel="log of difference", titlefont=font(10,), title="Logarithm of the norm of difference of the orbidts");
savefig(diff_plot, "lorenz-figs/fig10.png");