📂力学系

カオス遷移

定義 ¹

ある一定の期間にわたってキャリック的に動いていた軌道が突然キャリック領域を逸脱することを 過渡的混沌^{transient chaos}という。過渡的混沌が持続する期間を 過渡的寿命^{transient lifetime}という。

説明 ²

食物連鎖システム

食物連鎖システム: $$ \begin{align*} \dot{R} =& R \left( 1 - {\frac{ R }{ K }} \right) - {\frac{ x_{c} y_{c} C R }{ R + R_{0} }} \\ \dot{C} =& x_{c} C \left( {\frac{ y_{c} R }{ R + R_{0} }} - 1 \right) - {\frac{ x_{p} y_{p} P C }{ C + C_{0} }} \\ \dot{P} =& x_{p} P \left( {\frac{ y_{p} C }{ C + C_{0} }} - 1 \right) \end{align*} $$

この系でパラメータを $x_{c} = 0.4$, $y_{c} = 2.009$, $x_{p} = 0.08$, $y_{p} = 2.876$, $R_{0} = 0.16129$, $C_{0} = 0.5$ と置き $K = 1.0001$ としたとき、ある初期条件に対して次のような時間発展が見られる。キャリック的に振る舞っているような $P$ の個体数はある時点 $t = \tau$ を境に急激に減少するが、この時点 $\tau$ がすなわち過渡的寿命になる。

このような過渡的混沌は初期条件によって変わる。初期条件を変えて10万回のシミュレーションを行うと、その過渡的寿命が指数分布に従うことが分かる。

ロジスティック写像

別の例としてロジスティックファミリーを考えると、$g_{a} = ax (1-x)$ で作られる系は パラメータ $a$ によって変化する様子 を示し、$a=4$ のときカオス軌道を持つ ことが確認できる。すると次の問いはすなわち「$a>4$ のときはどうなるか」である。

まず $a>4$ ならば $g_{a}(x) = ax(1-x)$ であるから、その最大値は $1$ より大きく、$g_{a}$ は $x_{n} \notin [-1,1]$ を負に写す。また二次関数なので、一度でもある $\tau \in \mathbb{N}$ に対して $x_{ \tau } > 1$ になればその後は猛烈な速度で発散する。カオス軌道はその前は有界軌道でなければならないが発散したということはパラメータ $a$ の変化によりカオス性を失ったことを意味し、上の定義に従ってカオス遷移と呼べる。

この現象は $a>4$ の値が大きいほど速く起きるだろう。すると $x_{ \tau } > 1$ になるような $\tau$ は $a$ とどのような関係があるか推測できる。これを実際にシミュレーションで確認するのはさほど難しくなく、コードを実行するだけで次のプロットが描かれる。

コードでは $a$ の値をごくわずかずつ増やしながらランダムな初期値 $x_{0}$ に対して $x_{n} >1$ を満たす $n = \tau$ を見つけ、それを $N=1000$ 回繰り返して $\tau$ たちの平均 $\overline{\tau}$ を計算する。実際に描かれた図を見ると $a = 4 + \varepsilon$ に対して $\log \varepsilon \sim \log \left( \overline{ \tau } \right)$ のような線形関係が現れることが確認できる。

このようなカオス遷移はロジスティック写像の系だけに現れるのではなく、多くの系で見られる現象だ。

コード

Julia

以下は食物連鎖システムで過渡的混沌を再現するためのJuliaコードだ。

function RK4(f::Function, v::AbstractVector, h=1e-2)
    V1 = f(v)
    V2 = f(v + (h/2)*V1)
    V3 = f(v + (h/2)*V2)
    V4 = f(v + h*V3)
    return v + (h/6)*(V1 + 2V2 + 2V3 + V4), V1
end

function factory_foodchain2(K::Number; ic = [0., 0.4rand() + 0.6, 0.4rand() + 0.15, 0.5rand() + 0.3], tspan = 0:1e-2:10)
     xc,    yc,   xp,    yp,      R0,  C0 = (
    0.4, 2.009, 0.08, 2.876, 0.16129, 0.5)
    function sys(v::AbstractVector)
        _,R,C,P = v

        dt = 1
        dR = R*(1 - frac(R,K)) - xc*yc*frac(C*R, R + R0)
        dC = xc*C*(frac(yc*R, R + R0) - 1) - xp*yp*frac(P*C, C + C0)
        dP = xp*P*(frac(yp*C, C + C0) - 1)
        return [dt, dR, dC, dP]
    end

    len_t_ = length(tspan); h = tspan.step.hi
    v = ic; DIM = length(v)
    traj = zeros(2+len_t_, 2DIM); traj[1, 1:DIM] = v

    tk = 0
    while tk ≤ len_t_
        t,_,_,_ = v
        v, dv = RK4(sys, v, h)

        if t ≥ first(tspan)
            tk += 1
            traj[tk+1,         1:DIM ] =  v
            traj[tk  , DIM .+ (1:DIM)] = dv
        end
    end
    return traj[1:(end-2), :]
end
factory_foodchain2(T::Type, args...; kargs...) =
DataFrame(factory_foodchain2(args...; kargs...), ["t", "R", "C", "P", "dt", "dR", "dC", "dP"])[:, Not(:dt)]

hrzn = []
vrtl = []
@time for k in 0.9:0.0001:1
    for _ in 1:500
        traj = factory_foodchain2(DataFrame, k, tspan = 0:1e-1:5000)[40000:end,:]
        bits = arglmin(traj.P)
        if maximum(traj.P) > 0.55
            push!(hrzn, fill(k, length(bits)))
            push!(vrtl, traj.P[bits])
            break
        else
            println("k: $k should be tried again")
        end            
    end
end
using LaTeXStrings
scatter([hrzn...;], [vrtl...;], ylims = [0.55, 0.8], yticks = 0.55:0.05:0.8, xticks = 0.9:0.02:1,
    ms = 1, msw = 0, color = :black, alpha = 0.5, legend = :none, xlabel = L"K", ylabel = L"\min P")
png("bifurcaiton.png")

tau_ = zeros(Int64, 100000)
@showprogress @threads for k in eachindex(tau_)
    traj = factory_foodchain2(DataFrame, 0.99976 + 2e-4, ic = rand(4), tspan = 0:1e-0:5000)
    tau = findlast(diff(cumsum(traj.P .< 0.55)) .== 0)
    tau_[k] = !isnothing(tau) ? tau : -1
end
freq = [count(x .≤ tau_ .< x + 250) for x in 250:250:4750]
scatter(250:250:4750, freq ./ sum(freq), yscale = :log10, xticks = 0:1000:5000, xlims = [0, 5000], yticks = [1e-1, 1e-2], color = :white, legend = :none, xlabel = "transient lifetime", ylabel = "frequency")

R

以下はロジスティック写像で過渡的混沌を再現するためのRコードだ。

ga<-function(a,x,k=1) {
  tmp <- a*x *(1 - x)
  if(k==1) {return( tmp ) }
  else {return(ga(a,tmp,k-1))}
}
 
N<-1000
y<-numeric(0)
A<-4+10^(-(80:120)/20)
for(a in A){
  taubar<-0
  for(i in 1:N){
    tau<-0
    x<-runif(1)
    while(x>0){
      x<-ga(a,x)
      tau=tau+1
    }
    taubar<-taubar+tau
  }
  y<-c(y,taubar/N)
}
 
taubar<-y
a<-A
 
win.graph(); plot(log(a-4),log(taubar),type='l',main='a>4 일 때 로지스틱 맵의 Chaotic Transient')