📂力学系

ブアリ・アトラクター

モデル¹

次はブアリ・アトラクター^{Bouali attractor}の支配方程式である。

$$ \begin{align*} \dot{x} =& x ( a - y ) + \alpha z \\ \dot{y} =& - y ( b - x^2 ) \\ \dot{z} =& - x ( c - s z ) - \beta x \end{align*} $$

パラメータ

• $a = 4$
• $\alpha = 0.3$
• $b = 1$
• $c = 1.5$
• $\beta = 0.05$
• $s = 1$

説明

$$ \begin{align*} \dot{x} =& x ( a - y ) \\ \dot{y} =& - y ( b - x^2 ) \end{align*} $$

ブアリ・アトラクター^{Bouali attractor}は、上のようなロトカ・ヴォルテラモデル型のシステムに$z$軸を追加してカオス化^{chaotification}したモデルである。$\dot{z}$の形が必ずしもあのようでなければならないわけではなく、自然な振動^oscillationを持って被食者である$x$にフィードバックを与えさえすればよい。

実際、ブアリ方程式といえば、もう一つの類似した形として次が知られており、典型的なパラメータは$\alpha = 3$、$\beta = 2.2$、$\gamma = 1$、$\mu = 0.001$である²。

$$ \begin{align*} \dot{x} =& a x (1 - y) - \beta z \\ \dot{y} =& - \gamma y (x - 1) \\ \dot{z} =& \mu x \end{align*} $$

アトラクター

$s = 1$のときのストレンジ・アトラクターである。ローレンツ・アトラクターと比べると、翼の下へ伸びていく何かが生じたような印象を受ける。

分岐

$s$を0.5から1.2まで変化させながら描いた分岐ダイアグラムである。最高次の項が3次だからか、相沢アトラクターに劣らず奇怪な形を見せている。

コード

次はブアリ・アトラクターをシミュレーションできるJuliaコードである。

using DataFrames, DifferentialEquations, OrdinaryDiffEqLowOrderRK

function factory_bouali(s::Number; ic = [1.0, 1.0, 0.1], saveat = 0:1e-2:100)
      a,   α, b,   c,    β = (
      4, 0.3, 1, 1.5, 0.05)
    function sys(du, u, p, t)
        x, y, z = u
        s = p[1]

        du[1] = x*(a - y) + α*z
        du[2] = -y*(b - x^2)
        du[3] = -x*(c - s*z) - β*z
        return du
    end
    sol = solve(ODEProblem(sys, ic, (0, last(saveat)), [s]), RK4(), dt = saveat.step.hi, adaptive=false, maxiters = Inf)
    dsol = [sol.k[2][1] stack([solk[2] for solk in sol.k[2:end]])]
    matrix = [sol.t sol[:, :]' dsol']
    return matrix[sol.t .≥ first(saveat), :][1:end-1, :]
end
factory_bouali(T::Type, args...; kargs...) =
DataFrame(factory_bouali(args...; kargs...), ["t", "x", "y", "z", "dx", "dy", "dz"])

sol = factory_bouali(DataFrame, 1, ic = [2,1,0], saveat = 4500:1e-2:5000)
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.5, 1.2; p0, p1 = 0.7, 0.8;
p_ = range(pm, pM, length = 1001)
bfcn = callbfcn("G:/BF/bouali/bfcnA.jld2")
@showprogress @threads for k in eachindex(p_)
    sol = factory_bouali(DataFrame, p_[k], saveat = 2500:1e-2:3000)
    bfcn[p_[k]] = sol.x[arglmax(sol.x)]
end
scatter(dict2bifurcation(bfcn)..., xticks = [pm, p0, p1, pM], ms = .5, ma = .5, msw = 0, color = :black); png("temp1")