📂力学系

アイザワアトラクター

モデル

$$ \begin{align*} \dot{x} =& (z - b)x - d y \\ \dot{y} =& d x + (z - b)y \\ \dot{z} =& c + a z - \frac{z^{3}}{3} - (x^{2} + y^{2})(1 + e z) + f z x^{3} \end{align*} $$

パラメータ

• $a = 0.95$
• $b = 0.7$
• $c = 3.5$
• $d = 0.25$
• $e = 0.1$
• $f = 0.1$

説明

アイザワアトラクター^{Aizawa attractor}は右辺に最大4次項までを含み複雑なカオスを見せる力学系であり、時系列予測手法などでベンチマークとして登場することもある。

しかし筆者はこのシステムの起源と由来を調べていくうちに、アイザワ本人の論文が存在せず、かろうじてアイザワアトラクターが登場する文献においても方程式の表現がまちまちであることを知った^{1 2 3 4 5}。すでに名声を得て引用されているアトラクターにしてはかなり奇妙なことである。本文では常微分方程式をそれらの中で最も普遍的な形態と思われる表現にまとめて残した。

アトラクター

ストレンジアトラクターの中でもかなり奇妙な形をしている。基本的に球の形態を取りつつも$z$軸に沿って回転しながら上下する構造を持っている。

バイファケーション

少なくともこのポストを書いている今に至るまで、筆者はアイザワアトラクターに関連するバイファケーションに言及した研究結果を見たことがない。上のバイファケーション図は$b$を$0.6$から$0.9$まで変化させながら、$z$の最大値を記録したものである。見ての通り、普通なら見えるはずのルート・トゥ・カオス^{route to chaos}などが見られず、非常に奇怪な形態を持っている。私は今はアイザワアトラクターを研究する時間がないが、いつかはこのシステムについて深く掘り下げてみたい。

コード

以下はアイザワアトラクターをシミュレーションできるJuliaコードである。

using DataFrames, DifferentialEquations, OrdinaryDiffEqLowOrderRK

function factory_aizawa(b::Number; ic = [0.1, 0.1, 0.1], saveat = 0:1e-2:100)
      a,  c,   d,   e,  f = (
    .95, .6, 3.5, .25, .1)
    function sys(du, u, p, t)
        x, y, z = u
        b = p[1]

        du[1] = (z - b)*x - d*y
        du[2] = d*x + (z - b)*y
        du[3] = c + a*z - (z^3)/3 - (x^2 + y^2)*(1 + e*z) + f*z*x^3
        return du
    end
    sol = solve(ODEProblem(sys, ic, (0, last(saveat)), [b]), RK4(), dt = saveat.step.hi, adaptive=false, maxiters = 1e+7)
    matrix = Matrix([sol.t'; sol[:, :]; stack([sys(zeros(3), u, [b], 0) for u in sol.u])]')
    return matrix[sol.t .≥ first(saveat), :][1:end-1, :]
end
factory_aizawa(T::Type, args...; kargs...) =
DataFrame(factory_aizawa(args...; kargs...), ["t", "x", "y", "z", "dx", "dy", "dz"])

sol = factory_aizawa(DataFrame, 0.7)
plot(sol.x, sol.y, sol.z, alpha = .5, color = :black, xlabel = L"x", ylabel = L"y", zlabel = L"z", size = [400, 400])

pm, pM = 0.6, 0.9; p0, p1 = 0.72, 0.74;
p_ = range(pm, pM, length = 1001)
bfcn = callbfcn("G:/BF/aizawa/bfcnA.jld2")
@showprogress @threads for k in eachindex(p_)
    sol = factory_aizawa(DataFrame, p_[k], saveat = 0:1e-3:3000)
    z_ = sol.z[sol.t .≥ 2000]
    bfcn[p_[k]] = z_[arglmax(z_)]
end