📂동역학

과도적 혼돈

정의 ¹

일정한 기간동안 캐어릭하게 움직이던 트래젝터리가 갑자기 캐어릭 영역을 벗어나는 것을 과도적 혼돈^{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)$ 과 같은 선형관계가 나타남을 확인할 수 있다.

이와 같은 카오틱 트랜지션은 로지스틱 맵의 시스템에서만 나타나는 것이 아니라 많은 시스템에서 발견되는 현상이다.

코드

줄리아

다음은 먹이사슬 시스템에서 과도적 혼돈을 재현하기 위한 줄리아 코드다.

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')