Hénon写像

記事の内容


今回は, Hénon写像で図を描いて遊びます. 時間ができたら, 不動点とか分岐の理論的な部分もやりたいな...

Hénon写像の基本

定義

今回の主役はHénon写像です. Hénon写像は2次元の離散力学系で, 定義は以下の通りです.

【定義: Hénon写像】

\(a,b\in \mathbb{R}\)に対して, Hénon写像\(F:\mathbb{R}^2\to \mathbb{R}^2\)を次式で定める:

\begin{equation} F\left( \begin{bmatrix} x_n \\ y_n \end{bmatrix} \right) = \begin{bmatrix} y_n+1-ax_n^2 \\ bx_n \end{bmatrix} \end{equation}

式を見てもよく分からないので, 図を描いてみます. aとbはしばしば, それぞれ1.4と0.3として選択されます.

軌道

適当な初期値から始めて, その後の時間発展を追います. つまり軌道をプロットします. 下図は, 原点を初期値とし, a=1.4, b=0.3とした場合の軌道です.

【コード3の実行結果】

放物線のようにぐにゃっと曲がっています. Hénon写像は次のように分解することができます.

\begin{equation} F = f_3 \circ f_2 \circ f_1 \end{equation}

ここで, \begin{align} f_1\left( \begin{bmatrix} x_n \\ y_n \end{bmatrix} \right) &= \begin{bmatrix} x_n\\ 1-ax_n^2+y_n \end{bmatrix} \\ f_2\left( \begin{bmatrix} x_n \\ y_n \end{bmatrix} \right) &= \begin{bmatrix} b x_n\\ y_n \end{bmatrix} \\ f_3\left( \begin{bmatrix} x_n \\ y_n \end{bmatrix} \right) &= \begin{bmatrix} y_n\\ x_n \end{bmatrix} \end{align} です. 引数のベクトルは, まずy軸方向に放物線になるよう曲げられます(\(f_1\)による変形). 次に, x軸方向にb倍に圧縮されます(\(f_2\)による変形). 最後に, y=xに関して対称に写されます(\(f_3\)による変形).

分岐

分岐の様子を観察します. パラメータが2つあるので, 少し工夫が必要です.

1パラメータの分岐図

まず, 一方のパラメータを固定して, 他方のパラメータに関する分岐を調べます. 最初に, bを0.3に固定して, aだけを変化させます. 分岐図の描き方はロジスティック写像の時と同じです. まず, 候補となるパラメータ値aを一つ固定し, 初期値を原点とします. Hénon写像を1000回ほど作用させ, 軌道がある程度落ち着くまで計算します. 次に, その軌道の最後の値から出発して, 追加で1000回分計算します. その1000個の値を, 予め固定したaに対して縦軸方向にプロットします. 下図はこうして得られた図です. 各aに対して, 縦方向に1000個分の点が描かれています.

【コード4の実行結果】

ロジスティック写像の分岐図に似ています. 周期倍分岐が起きている様子も観察できます.

次に, aを1.4に固定して, bだけを変化させます. 描き方は先ほどと同様です.

【コード5の実行結果】

こちらも, ロジスティック写像の分岐図に大変似ています.

2パラメータの分岐図

前節では, パラメータのうち一方を固定し, 他方を変化させていました. 今度は, どちらの値も変化させます. ここでは, 文献[2]の方法に従って, 周期軌道を表現します.

まず, パラメータ値a,bを固定します. そして, 初期値を原点とします. Hénon写像を1000回ほど作用させ, 軌道がある程度落ち着くまで計算します. 次に, その軌道の最後の値(以下, "最後の点")から出発して, Hénon写像を最大\(N\)回ほど作用させ, 軌道の周期性を判定します. ただし, 写像を1回作用させるごとに判定します. 写像を1回作用させた"予測値"と, "最後の点"との距離を\(d\)とします. \(d\)が十分小さい値\(\varepsilon\)より小さければ, "最後の点"と同じ点を通過したとみなし, 周期点と判定します. そのときの反復回数を周期とします. \(d\)が十分大きい値\(R\)より大きければ発散したとみなし, 計算はストップします. 周期性はないと判定します. どちらでもなければ, 計算を続けます. こうして, 最大\(N\)周期までを発見することができます. このようにして周期を判定して描いた分岐図が下図です. \(\varepsilon=,10^{-6} R=10^6, N=32\)としました. 色が薄いところほど周期は小さいく, 色が濃くなるほど周期は大きくなります. 黒い部分は発散しています.

【コード6の実行結果】

コード

インポートと関数の定義

【Juliaコード1; インポート】
using Plots
using LinearAlgebra
pyplot()
【Juliaコード2; 関数の定義】
#Hénon写像
Hénon_map(x::Float64, y::Float64, a::Float64, b::Float64) = y+1-a*x^2, b*x

#Hénon写像の軌道
function Hénon_orbit(a::Float64, b::Float64, x₀::Float64, y₀::Float64, n_iter::Int64)
    xs = zeros(n_iter)
    ys = zeros(n_iter)
    xs[1] = x₀
    ys[1] = y₀
    for i in 2:n_iter
        xs[i], ys[i] = Hénon_map(xs[i-1], ys[i-1], a, b)
    end
    xs, ys
end

#周期Tを計算する
function compT(a::Float64, b::Float64, x::Float64, y::Float64, N::Int64, ϵ::Float64, R::Float64)
    #周期
    T = 0

    #新たな点
    x_new = x
    y_new = y

    #N回計算する
    for n in 1:N
        #Hénon写像で写す
        x_new, y_new = Hénon_map(x_new, y_new, a, b)

        #最後の点との距離
        d = sqrt((x-x_new)^2+(y-y_new)^2)

        #発散する場合
        if d>R
            T = -1
            return T
        #近傍を通る場合
        elseif d<ϵ
            return T
        #周期点でもなく, 発散もしない場合
        else
            T += 1
        end

        #周期点でない場合
        if T>N
            T = -1
        end
    end
    return T
end

軌道

【Juliaコード3; 軌道の計算】
#a=1.4, b=0.3
a = 1.4
b = 0.3
x₀ = 0.0
y₀ = 0.0
xs, ys = Hénon_orbit(a,b,x₀,y₀,10000)

#写像の軌道をプロット
fig1 = plot(title="Hénon Map: a=$(a), b=$(b)", xlabel="xn", ylabel="yn")
plot!(xs, ys, st=:scatter, markersize=1, label=false)
savefig(fig1, "figs-Henon_map/fig1.png")

分岐図

【Juliaコード4; 1パラメータの分岐図】
#b=0.3
b = 0.3

#初期値
x₀ = 0.0
y₀ = 0.0

#反復回数, aの候補, 候補数
n_iter = 1000
as = range(0.0, 1.5, length=500)

#保存用
xs = zeros(n_iter)
ys = zeros(n_iter)

#軌道の落ち着き先の計算し, プロット
fig2 = plot(title="Bifurcation: Hénon Map", xlabel="a", ylabel="xn")
for a in as
    xs, ys = Hénon_orbit(a,b,x₀,y₀,n_iter)
    xs, ys = Hénon_orbit(a,b,xs[end],ys[end],n_iter)
    plot!([a], xs, st=:scatter, markersize=0.1, label=false)
end
savefig(fig2, "figs-Henon_map/fig2.png")
【Juliaコード5; 1パラメータの分岐図】
#a=1.4
a = 1.4

#初期値
x₀ = 0.0
y₀ = 0.0

#反復回数, aの候補, 候補数
n_iter = 1000
bs = range(-0.5, 0.3, length=500)

#保存用
xs = zeros(n_iter)
ys = zeros(n_iter)

#軌道の落ち着き先の計算し, プロット
fig3 = plot(title="Bifurcation: Hénon Map", xlabel="b", ylabel="xn")
for b in bs
    xs, ys = Hénon_orbit(a,b,x₀,y₀,n_iter)
    xs, ys = Hénon_orbit(a,b,xs[end],ys[end],n_iter)
    plot!([b], xs, st=:scatter, markersize=0.1, label=false)
end
savefig(fig3, "figs-Henon_map/fig3.png")
【Juliaコード6; 2パラメータの分岐図】
#N, ϵ, Rは固定する
N = 32
ϵ = 1e-6
R = 1e6

#aの候補
as = range(0.0, 2.0, length=600)
na = length(as)

#bの候補
bs = range(-1.0, 1.0, length=600)
nb = length(bs)

#xとyの保存用
n_iter = 1000
x = 0.0
y = 0.0

#周期の保存用
Ts = zeros(Int64, nb, na)

#aとbを変化させて計算する
for j in 1:nb
    b = bs[j]

    for i in 1:na
        a = as[i]

        #初期化
        x = 0.0
        y = 0.0

        #初期値から初めて軌道を計算する
        for itr in 1:n_iter
           x,y = Hénon_map(x,y,a,b)
        end

        #周期を計算する
        T = compT(a, b, x, y, N, ϵ, R)
        Ts[j,i] = T
    end
end
参考文献
      [1]S.H.Strogatz, 田中久陽, 中尾裕也, 千葉逸人, ストロガッツ 非線形ダイナミクスとカオス, 丸善出版, 2015
      [2]小室 元政, 基礎からの力学系―分岐解析からカオス的遍歴へ, サイエンス社, 2005