ロジスティック写像まとめ

記事の内容


今回はロジスティック写像を解析します. 機会があれば, エントロピー, クライシスなども追加します. とりあえず今回は, カオス的な振る舞いや, 周期的な振る舞いなど, 図を描いてわかる範囲で記事にします.

ロジスティック写像の概要

定義

ここでは, 写像という語を, 普通の意味とは異なった意味で用います.

【定義: 写像/離散時間力学系】

連続写像\(f:\mathbb{R}\to \mathbb{R}\)に対して, 次の差分方程式を離散時間力学系, または写像という.

\begin{equation} x_{n+1} = f(x_n) \quad (n=0,1,\cdots) \end{equation}

また, \(\{x\}_{n=0}^{\infty}\)を軌道という.

それでは, メインのロジスティック写像を定義しておきましょう. ロジスティック方程式と同じ形をしています.

【定義: ロジスティック写像】

ロジスティック写像を次のように定義する. ただし, \(r>0\)とする.

\begin{equation} x_{n+1} = rx_n(1-x_n) \end{equation}

以下, 次の記号を用います.

\begin{equation} f(x) = rx(1-x) \end{equation}

また, \(r\)のことをパラメータと呼びます.

解析の基本的な道具

写像の解析に際して, 次の解析ツールは基本的です.

  • 定性的な解析のための道具
    • 時系列プロット
    • クモの巣プロット
    • 軌道図
  • 定量的な解析のための道具
    • 線形安定性解析
    • Lyapunov指数
    • エントロピー

時系列プロットは, \(f\)を次々に作用させていった値を順にプロットしたものです. 直感的にも分かりやすく, 規則性の発見に役立ちます.

クモの巣プロットは, 次の順にプロットしていきます.

【アルゴリズム: クモの巣プロット】

\((x_0, 0)\)から\((x_0, x_1)\)に垂直な線分を引く.
for \(k=1,\cdots,K-1\)
\((x_{k-1}, x_k)\)から\((x_k, x_k)\)まで水平線を引く.
\((x_k, x_k)\)から\((x_k, x_{k+1})\)まで垂直線を引く.

クモの巣プロットを見れば, 不動点や周期点の存在, それらの安定性などが分かります.

軌道図は, 横軸にパラメータ値, 縦軸に十分な回数計算した軌道の値をプロットしたものです. 各パラメータごとに, 軌道の落ち着き先が分かります. ロジスティック写像では, 非常に特徴的な図になります.

定性的な解析

下図は, ロジスティック写像の有名な軌道図です. 横軸は, パラメータ\(r\), 縦軸は軌道が落ち着いた先の値で, 2000点からなる軌道の最後の100個の値を表します. 初期値は0.01です. 以下, 初期値は全て0.01です.

【コード4の実行結果】

上の軌道図にはいくつか特徴があります. 例えば, 左側の\(3\leq r \leq 3.5\)付近では, 軌道が2つの値に分かれています. この部分では, 2周期で値が写り合います. このような点を2周期点といいます. この用語をもう少し厳密に定義します.

【定義: p周期点】

写像\(x_{n+1}=f(x_n)\)に対して, 次を満たす点\(x^*\)をこの写像の\(p\)周期点という.

\begin{equation} f^p(x^*) = x^* ,\quad f^i(x^*) = x^*\ (i=1,\cdots,p-1) \end{equation}

\(p\)周期点を含む\(p\)周期軌道という.

2周期軌道

分岐図でまず目に付く2周期軌道を観察していきましょう. 下図は, 分岐図の2周期軌道部分を取り出したものです. ちょうど\(r=3.0\)で二股に分かれています. ちょうど\(3.0\)であることは, 解析的に示すことができます.

【コード5の実行結果】

次に, 時系列をプロットしておきます. 早い段階で, 2周期軌道に落ち着きます. 以下, \(r=3.2\)です.

【コード6の実行結果】

下図は, 蜘蛛の巣プロットです. 200反復分計算しています. 時系列プロットに対応して, こちらも早い段階で2周期の軌道に落ち着いています. 線が何度も重なっている部分が周期軌道です.

【コード7の実行結果】

下図は, クモの巣プロットのアニメーションです. こちらは100反復分計算しています. 2周期の軌道に巻きついていく様子が観察できます.

【コード8の実行結果】

3周期軌道

次に3周期軌道を計算します. 下図は, 軌道図の一部を切り出したものです. 縦に3個に分かれているので, これらの値しか取らない3周期軌道が存在すると考えられます.

【コード9の実行結果】

時系列で見ると, 確かに3周期軌道です. 以下, \(r=3.83\)です.

【コード10の実行結果】

下図は, 200反復計算したクモの巣プロットです. 少し分かりにくいですが, 黒線が重なった周期軌道が見えます.

【コード11の実行結果】

下図は, 100反復計算したクモの巣プロットです. こちらは周期軌道がはっきりと現れています. 上図より100反復分少ないので, あと100反復のうちに, 一度周期性が崩れることが示唆されます.

【コード12の実行結果】

カオス的な挙動

最後に, 周期軌道ではなく, 軌道図で黒く塗りつぶされた部分を見てみます.

【コード13の実行結果】

時系列で見ると, 周期性は観察されません. 以下, \(r=3.752\)です.

【コード14の実行結果】

蜘蛛の巣プロットも複雑です. 200反復計算しています. デタラメな値をとっているように見えます.

【コード15の実行結果】

下図は100反復分計算したアニメーションです.

【コード16の実行結果】

定量的な解析

ここまでは主にグラフを見ながら話を進めてきました. ここではもう少し式や具体的な数値を交えて話を進めていきます.

線形安定性解析

線形安定性解析は, 微分係数から不動点の安定性を判定する方法です.

まず, \(x^*\)を不動点とします. 軌道の不動点からのずれを\(\xi_n\)で表します. すなわち, 次のように定めます.

\begin{equation} \xi_n = x_n - x^* \end{equation}

イメージ的には, 不動点近傍の軌道を引き寄せる, 吸引的な不動点を, 安定な不動点といいます. 逆に, 不動点近傍の軌道を遠ざける, 反発的な不動点を, 不安定な不動点といいます. したがって, \(\xi_n\)が減衰するならば安定, 増幅するならば不安定だと判定できそうです.

それでは, 減衰・増幅の判定方法を考えます. 不動点の近傍で, 次のようにTaylor展開します.

\begin{align} x_{n+1} &= f(x_n) \\ &= f(x^* + \xi_n) \\ &= f(x^*) + f^\prime(x^*) \xi_n + O(\xi_n^2)\\ &= x^* + f^\prime(x^*) \xi_n + O(\xi_n^2) \end{align}

したがって, 2次以上の項を無視して, 次式が成り立ちます.

\begin{equation} \xi_{n+1} = f^\prime(x^*) \xi_n \end{equation}

以上より, 次式が成り立ちます.

\begin{equation} \xi_{n} = (f^\prime(x^*) )^n\xi_0 \end{equation}

この式から, 次のような挙動が判明します.

  • \(f^\prime(x^*) < 1\)ならば, \(\xi_{n}\to 0 \ (n\to \infty)\)である. つまり, 不動点は安定である.
  • \(f^\prime(x^*) > 1\)ならば, \(\xi_{n}\to \infty \ (n\to \infty)\)である. つまり, 不動点は不安定である.

微分係数がぴったり1の場合には, 非線形項の影響を加味する必要があります.

2周期軌道

ここでは, 2周期軌道がどのようにして発生するかを解析します. 2周期軌道は, \(f^2\)が対角線\(y=x\)で交わるところに発生します. つまり, \(f^2\)の不動点です. しかし, \(f^2\)の不動点には\(f\)の不動点が含まれます. \(f\)の不動点を除外する必要があります. 以下, \(f\)の2周期点を探します. 手順としては, 「\(f^2\)の不動点を列挙する」, そして, 「\(f\)の不動点を除外する」, の2段階です. なお, 以下\(r>1\)としておきます.

その前に, 2周期点と2周期軌道が何個存在するかを数えます. まず, \(f\)には2個の不動点が存在します. また, これを2回作用させた\(f^2\)には, \(2^2=4\)個の不動点が存在します. 一方, \(f\)の不動点が2個存在するので, 4個中2個は\(f\)の不動点です. したがって, \(f\)の2周期点は残りの2個です. この2個の不動点が\(f\)により互いに写りあうので, \(f\)の2周期軌道は1個です. 以上より, パラメータ\(r\)を変化させるとどこかで2周期軌道が1つだけ発生します.

2周期点の値を解析的に求めるために, 次式を解きます.

\begin{equation} f^2(x) = x \end{equation}

この方程式の解は最大で4個存在します. しかし, 上で説明したように, そのうち2個は\(f\)の不動点です. \(f\)の不動点が \[ x = 0, 1-\frac{1}{r} \]なので, 以下の形になります.

\begin{equation} x\left\{ x-\left( 1-\frac{1}{r}\right)\right\} \left\{ r^3x^2 -(r+1)r^2x + (r^2+r)\right\} = 0 \end{equation}

結局, \(f\)の2周期点を求めるには, 後半の2次方程式を解けば良いです. \(f\)の2周期点は次の通りです.

\begin{equation} x = \frac{(r+1) \pm \sqrt{(r+1)(r-3)}}{2r} \end{equation}

この2周期点が存在するためには, \(r>3\)が必要です. 軌道図において\(r=3\)で軌道が二股に分かれるのはこのためです. これを別の角度から眺めてみます. 下図は, パラメータ\(r\)を大きくしたときの\(f^2\)の変化です. \(f^2\)の不動点が\(r=3\)で生じていることが分かります. この不動点の安定性はどうなっているでしょうか. \(r< 3\)では, \(f\)の右側の不動点における接線の傾きは1より小さいので, 安定です. これが\(r>3\)になると, \(f\)の右側の不動点における接線の傾きが1より大きくなり, 不安定になります. 一方, \(r=3\)を跨いで新たに生じた2つの不動点は, 接線の傾きから, 安定です. まとめると, 「パラメータ\(r\)を大きくすると, 既に存在した不動点が不安定化し, その周りに安定な不動点が生じる」ことが分かります. この現象を周期倍分岐といいます.

【コード18の実行結果】

3周期軌道

2周期軌道と同じ要領で考えます. まず, \(f^3\)の不動点は\(2^3=8\)個存在します. そのうち, \(f\)の不動点が2個なので, \(f\)の3周期点は残りの6個です. この6個で互いに3周期で写りあうので, 3周期軌道は2個です.

3周期軌道の発生の仕方は, 2周期軌道の時とは少し違います. 今回は不動点が突然出現します. グラフが対角線に触れた次の瞬間, その両側に安定な不動点と不安定な不動点が同時に出現します. このように, 「パラメータを大きくすると, 安定な不動点と不安定な不動点が同時に発生する」分岐をサドル・ノード分岐といいます.

【コード19の実行結果】

Lyapunov指数

写像の複雑さを定量評価する指標を用意します. \(n\)を1つ固定します. \(f\)を\(n\)回施したとき, 初期値\(x_0\)からの微小のズレ\(\delta\)がどの程度増幅・減衰するかを評価します. Taylor展開により, 次式が成り立ちます.

\begin{equation} f^n(x_0+\delta) - f^n(x_0) = \left. \frac{d}{dx}f^n(x)\right|_{x=x_0} \delta + O(\delta^2) \end{equation}

ここで, 右辺の微分係数は, 次のように表されます.

\begin{equation} \left. \frac{d}{dx}f^n(x)\right|_{x=x_0} = \prod_{i=0}^{n-1} f^\prime (x_i) \end{equation}

したがって, 微小のズレの減衰・増幅率として, 次の指標が意味を持ちます.

【定義: Lyapunov指数】

写像の軌道\(\{x_i\}_{i=0}^{\infty}\)に対して, Lyapunov指数を次式で定める.

\begin{equation} \lim_{n\to \infty} \frac{1}{n}\sum_{i=0}^{n-1} \log |f^\prime (x_i)| \end{equation}

Lyapunov指数が正ならば, 初期値の微小なズレは増幅されたと考えられます. Lyapunov指数が負ならば, 初期値の微小なズレは減衰したと考えられます. 下図では, 安定な周期軌道が生じている部分ではLyapunov指数が負になっています.

【コード20の実行結果】

ソースコード

共通のコード

【Juliaコード1; インポート】
using Plots
using Statistics
using Formatting
pyplot()
【Juliaコード2; ロジスティック写像の定義】
#ロジスティック写像
logistic_map(x::Float64, r::Float64) = 0≤x≤1 ? r*x*(1-x) : NaN

#初期値
x₀=0.01

定性的な解析で用いたコード

【Juliaコード3;関数の定義】
#時系列プロット
function time_series_plot(r::Float64, x₀::Float64)
    #軌道の長さ, 軌道の保存
    ts_length = 80
    ts_array = zeros(ts_length)
    ts_array[1] = x₀
    ts_plot = plot(title="time series", xlabel="n", ylabel="x_n", legend=false)

    #軌道の計算
    for k in 2:ts_length
        ts_array[k] = logistic_map(ts_array[k-1], r)
    end

    #プロット
    plot!(1:ts_length, ts_array, marker=true, markerstrokealpha=0.5)
    return ts_plot
end

#蜘蛛の巣プロット作成用
function radar_chart(r::Float64, x₀::Float64, n::Int64)
    #変数の初期化と初期値の保存
    x = x₀
    y = 0.0
    radar_chart_array = zeros(2,2*n+1)
    radar_chart_array[:,1] = [x; y]

    #rを固定したロジスティック写像
    my_log_map(x) = logistic_map(x,r)

    for k in 1:n
        #縦に, グラフとぶつかるまで移動
        y = logistic_map(x,r)
        radar_chart_array[:,2*k] = [x; y]

        #横に, 対角線とぶつかるまで移動
        x = y
        radar_chart_array[:,2*k+1] = [x; y]
    end

    #対角線と連続版の写像プロット
    radar_chart_plot = plot(0.0:0.01:1.0, my_log_map, color=:black, legend=false, xlabel="xn", ylabel="x_n+1")
    plot!(0:1, 0:1, color=:black)
    xfrom = radar_chart_array[1,1]
    yfrom = radar_chart_array[2,1]
    for i in 1:n
        #縦線のプロット
        xto = radar_chart_array[1, 2*i]
        yto = radar_chart_array[2, 2*i]
        plot!([xto,], [yto,], st=:scatter, color=:black)
        plot!([xfrom, xto], [yfrom, yto], color=:black)

        #更新
        xfrom = xto
        yfrom = yto

        #横線のプロット
        xto = radar_chart_array[1, 2*i+1]
        yto = radar_chart_array[2, 2*i+1]
        plot!([xto,], [yto,], st=:scatter, color=:black)
        plot!([xfrom, xto], [yfrom, yto], color=:black)

        #更新
        xfrom = xto
        yfrom = yto
    end
    return radar_chart_plot
end

#蜘蛛の巣プロット作成用
function radar_chart_anim(r::Float64, x₀::Float64, n::Int64)
    #変数の初期化と初期値の保存
    x = x₀
    y = 0.0
    radar_chart_array = zeros(2,2*n+1)
    radar_chart_array[:,1] = [x; y]

    #rを固定したロジスティック写像
    my_log_map(x) = logistic_map(x,r)

    for k in 1:n
        #縦に, グラフとぶつかるまで移動
        y = logistic_map(x,r)
        radar_chart_array[:,2*k] = [x; y]

        #横に, 対角線とぶつかるまで移動
        x = y
        radar_chart_array[:,2*k+1] = [x; y]
    end

    #対角線と連続版の写像プロット
    radar_chart_plot = plot(0.0:0.01:1.0, my_log_map, color=:black, legend=false, xlabel="xn", ylabel="x_n+1")
    plot!(0:1, 0:1, color=:black)
    xfrom = radar_chart_array[1,1]
    yfrom = radar_chart_array[2,1]
    anim = @animate for i in 1:n
        #縦線のプロット
        xto = radar_chart_array[1, 2*i]
        yto = radar_chart_array[2, 2*i]
        plot!([xto,], [yto,], st=:scatter, color=:black)
        plot!([xfrom, xto], [yfrom, yto], color=:black)

        #更新
        xfrom = xto
        yfrom = yto

        #横線のプロット
        xto = radar_chart_array[1, 2*i+1]
        yto = radar_chart_array[2, 2*i+1]
        plot!([xto,], [yto,], st=:scatter, color=:black)
        plot!([xfrom, xto], [yfrom, yto], color=:black)

        #更新
        xfrom = xto
        yfrom = yto
    end
    return anim
end

#軌道図
function orbit(x₀::Float64, r_min::Float64, r_max::Float64)
    #軌道の保存用
    orbit_length = 1000
    orbit_array = zeros(orbit_length)
    orbit_array[1] = x₀

    rs = range(r_min, r_max, length=2000)
    orbit_plot = plot(title="Bifurcation", xlabel="r", ylabel="xn", legend=false, xlim=[r_min,r_max])
    for r in rs
        #1000回分計算
        for k in 2:orbit_length
            orbit_array[k] = logistic_map(orbit_array[k-1], r)
        end

        #初期値の保存
        orbit_array[1] = orbit_array[end]

        #追加で1000回計算
        for k in 2:orbit_length
            orbit_array[k] = logistic_map(orbit_array[k-1], r)
        end

        #プロット
        plot!([r,], orbit_array[end-Int(orbit_length/10):end], st=:scatter, markersize=0.1)
    end
    return orbit_plot
end
【Juliaコード4: 軌道図】
fig1 = orbit(x₀, 2.8, 4.0)
savefig(fig1, "figs-logistic-map/fig1.png")
【Juliaコード5: 軌道図】
fig2= orbit(x₀, 2.8, 3.5)
savefig(fig2, "figs-logistic-map/fig2.png")
【Juliaコード6: 時系列プロット】
fig3 = time_series_plot(3.2, x₀)
savefig(fig3, "figs-logistic-map/fig3.png")
【Juliaコード7: クモの巣プロット】
fig4 = radar_chart(3.2, x₀, 200)
savefig(fig4, "figs-logistic-map/fig4.png")
【Juliaコード8: クモの巣プロットのアニメーション】
anim1 = radar_chart_anim(3.2, x₀, 100)
gif(anim1, "figs-logistic-map/anim1.gif", fps=20)
【Juliaコード9: 軌道図】
fig5 = orbit(x₀, 3.82, 3.87)
savefig(fig5, "figs-logistic-map/fig5.png")
【Juliaコード10: 時系列プロット】
fig6 = time_series_plot(3.83, x₀)
savefig(fig6, "figs-logistic-map/fig6.png")
【Juliaコード11: クモの巣プロット】
fig7 = radar_chart(3.83, x₀, 200)
savefig(fig7, "figs-logistic-map/fig7.png")
【Juliaコード12: クモの巣プロットのアニメーション】
anim2 = radar_chart_anim(3.83, x₀, 100)
gif(anim2, "figs-logistic-map/anim2.gif", fps=20)
【Juliaコード13: 軌道図】
fig8 = orbit(x₀, 3.7, 3.8)
savefig(fig8, "figs-logistic-map/fig8.png")
【Juliaコード14: 時系列プロット】
fig9 = time_series_plot(3.752, x₀)
savefig(fig9, "figs-logistic-map/fig9.png")
【Juliaコード15: クモの巣プロット】
fig10 = radar_chart(3.752, x₀, 200)
savefig(fig10, "figs-logistic-map/fig10.png")
【Juliaコード16: クモの巣プロットのアニメーション】
anim3 = radar_chart_anim(3.752, x₀, 100)
gif(anim3, "figs-logistic-map/anim3.gif", fps=20)

定量的な解析で用いたコード

【Juliaコード17: Lyapunov指数のプロット用関数】
#蜘蛛の巣プロット作成用
function myLyapunov(x₀::Float64)
    #rの数
    rs_length = 2000

    #Lyapunov指数の保存用
    lyapunov_array = zeros(rs_length)

    #軌道の保存用
    orbit_length = 1000
    orbit_array = zeros(orbit_length)
    orbit_array[1] = x₀

    rs = range(2.8, 4.0, length=rs_length)
    for i in 1:rs_length
        r = rs[i]

        #1000回分計算
        for k in 2:orbit_length
            orbit_array[k] = logistic_map(orbit_array[k-1],r)
        end

        #初期値の保存
        orbit_array[1] = orbit_array[end]

        #追加で1000回計算
        for k in 2:orbit_length
            orbit_array[k] = logistic_map(orbit_array[k-1], r)
        end

        #Lyapunov指数の計算
        lyapunov_array[i] = mean(log.(abs.(r*(1 .- 2*orbit_array))))
    end
    lyapunov_plot = plot(
        rs,
        lyapunov_array,
        st=:scatter,
        markersize=1,
        title="Lyapunov exponent",
        xlabel="r",
        ylabel="Lyapunov exponent",
        legend=false
    )
    return lyapunov_plot
end
【Juliaコード18: 2周期点のアニメーション】
logmap2(x,r) = logistic_map(logistic_map(x,r),r)
rs = range(2.6, 3.4, length=100)
anim4 = @animate for r in rs
    plot(
        0:0.01:1,
        x->logmap2(x,r),
        title="r=$(format(r, precision=2))",
        xlabel="x",
        xlim=[0,1],
        ylim=[0,1],
        label="f^2"
    )
    plot!(0:1, 0:1, label="y=x")
end
gif(anim4, "figs-logistic-map/anim4.gif", fps=10)
【Juliaコード19: 3周期点のアニメーション】
logmap3(x,r) = logistic_map(logistic_map(logistic_map(x,r),r),r)
rs = range(3.81, 3.88, length=100)
anim5 = @animate for r in rs
    plot(
        0:0.005:1,
        x->logmap3(x,r),
        title="r=$(format(r, precision=3))",
        xlabel="x",
        xlim=[0,1],
        ylim=[0,1],
        label="f^3",
        legend=:bottomright
    )
    plot!(0:1, 0:1, label="y=x")
end
gif(anim5, "figs-logistic-map/anim5.gif", fps=10)
【Juliaコード20: Lyapunov指数のプロット】
fig11 = myLyapunov(x₀)
savefig(fig11, "figs-logistic-map/fig11.png")
参考文献
      [1]S.H.Strogatz, 田中久陽, 中尾裕也, 千葉逸人, ストロガッツ 非線形ダイナミクスとカオス, 丸善出版, 2015
      [2]小室 元政, 基礎からの力学系―分岐解析からカオス的遍歴へ, サイエンス社, 2005