力学系を観察しよう

記事の内容


今回はvan der Pol方程式を定性的に解析します.

2次元力学系: van der Pol方程式

van der Pol方程式

2次元力学系の例として, van der Pol方程式を考えます. 次式の通りです.

\[ \ddot{x} = (\mu-x^2)\dot{x} -x \]

2階微分方程式のままでは扱いにくいので, 次のように変形します. まず, xの1階微分をv表します. すると, van der Pol方程式は次のような微分方程式系になります.

\[ \frac{d}{dt}\begin{bmatrix} x\\ v \end{bmatrix} = \begin{bmatrix} v \\ -x^2v-x+\mu v \end{bmatrix} \]

力学系の解析では主に, 不動点やその周りでの軌道に関心があります. そこでこの方程式の不動点に注目します. van der Pol方程式の不動点は原点のみです.

線形化

非線形方程式を取り扱うのは一般に困難です. そこで不動点の周りで線形化します. 一般にはJacobi行列を計算しますが, 今回は式に形から簡単に求まります.

\[ \begin{bmatrix} v \\ -x^2v-x+\mu v \end{bmatrix} = \begin{bmatrix} 0 & 1\\ -1 & \mu \end{bmatrix} \begin{bmatrix} x\\ v \end{bmatrix} + \begin{bmatrix} 0\\ -x^2v \end{bmatrix} \] 不動点である原点付近では非線形項の寄与は小さいことを期待して, 次のように近似します. 固有値の値によっては, この近似は破綻します. \[ \frac{d}{dt}\begin{bmatrix} x\\ v \end{bmatrix} \simeq \begin{bmatrix} 0 & 1\\ -1 & \mu \end{bmatrix} \begin{bmatrix} x\\ v \end{bmatrix} \]

これは微分方程式の教科書に出てくる線形系です. 原点付近の軌道を求めるには, 適当に初期値を与えた初期値問題を解けば良いわけです. ただし, この系の解法には行列の指数関数などが出てきて, 複雑です. 一般に解は肩に固有値をのせた指数関数の一次結合です. よって, 固有値の正負, 相対的な大小が分かれば大まかな解析は可能です. 以下, 固有値を見ることで, van der Pol方程式を解析します.

固有値による分類と相図

線形系の形から, ベクトル場(相図, 相平面図)が描けます. また, 数値解法を利用すれば元の非線形系の軌道計算も簡単にできます. 以下, 固有値の値で分類して, 数値解, および相図を示します. 以下の図では, 赤丸が不動点, 青線が数値解, 黒矢印がベクトル場です. 数値解法として, 4段4次Runge-Kutta法を用いました. 反復回数は500回, 刻み幅は0.1です.

負の実根が2つのとき

このとき, μ<-2です. 固有値は負なので, 指数関数的に不動点に引き寄せられます. 横方向に横切るのが固有ベクトル方向です.

【コード3の実行結果】

負の実根が1つのとき

このとき, μ=-2です. 線形化した時の行列のランクが1つだけ落ちます. 原点付近では原点に巻きつくような形になります.

【コード4の実行結果】

実部が負の複素共役根のとき

このとき, -2<μ<0です. 原点に吸い寄せられるような軌道です.

【コード5の実行結果】

純虚数根のとき

このとき, μ=0です. 線形化したものと軌道がずれています. ベクトル場によれば不動点の周りを周期的に回転する軌道が予測されますが, 数値解は不動点に吸い寄せられています. これは線形系で無視した非線形項の影響です.

【コード6の実行結果】

実部が正の複素共役根のとき

このとき, 0<μ<2です. 軌道が巻きつくリミットサイクルが観察できます.

【コード7の実行結果】

正の実根が1つのとき

このとき, μ=2です. こちらもリミットサイクルが観察できます.

【コード8の実行結果】

正の実根が2つのとき

このとき, μ>2です. こちらもリミットサイクルが観察できます.

【コード9の実行結果】

ソースコード

共通のコード

【Juliaコード1; インポート】
using Plots
【Juliaコード2; 関数などの用意】
#右辺関数の線形化
linf(x,y,μ) = [0 1; -1 μ]*[x,y];

#右辺関数
f(x, μ) = [x[2]; -x[1]^2*x[2]-x[1]+μ*x[2]];

#グラフの端
e =4.0;

#RK法のパラメータ
n = 500; #反復回数
h = 0.1; #刻み幅
x₀_arr = [[-e/10, -e/10], [e/10, e/10], [-e; e], [-e,-e], [-1.0; -e], [0,e], [1.0; -e], [e,-e], [e; e], [-e, 0.0]]; #初期値

#4段4次RK法の1ステップ
function one_step_RK(x_old::Array{Float64,1}, h::Float64, f, μ::Float64)
    #kの更新
    k₁ = f(x_old, μ);
    k₂ = f(x_old.+h*k₁/2, μ);
    k₃ = f(x_old.+h*k₂/2, μ);
    k₄ = f(x_old.+h*k₃, μ);

    #qとvの更新
    x_new = x_old .+ h*(k₁.+2*k₂.+2*k₃.+k₄)/6;
    return x_new
end

#RK法によるシミュレーション
function myRK(n::Int64, x₀::Array{Float64,1}, f, μ::Float64)
    #変数の初期化
    xvec = zeros(length(x₀),n);
    xvec[:,1] = x₀;

    #反復計算
    for k in 1:n-1
        x_old = xvec[:,k];

        #qとvの更新
        x_new = one_step_RK(x_old, h, f, μ);

        #更新値の保存
        xvec[:,k+1] = x_new;
    end
return xvec
end

#プロット用のベクトルを作成する
function make_vectors(start::Union{Float64,Int64}, stop::Union{Float64,Int64}, step::Union{Float64,Int64}, f)
    #xとyの範囲
    xrange = range(start, stop, step=step);
    yrange = range(start, stop, step=step);

    #始点
    xstart = [x for x in xrange, y in yrange]';
    ystart = [y for y in yrange, x in xrange];

    #始点の配列
    xstart_array = xstart[:];
    ystart_array = ystart[:];

    #方向ベクトルのタプル
    vecfield = f.(xstart_array, ystart_array);
    xdir = [vecfield[i][1] for i in 1:length(vecfield)];
    ydir = [vecfield[i][2] for i in 1:length(vecfield)];
    quiver_tupple = (xdir, ydir)
    return xstart_array, ystart_array, quiver_tupple
end

相図の描画

【Juliaコード3; 負の実根が2つのとき】
#パラメータ値
μ = -3.0;

#ベクトル
xstart_array, ystart_array, quiver_tupple = make_vectors(-e, e, 0.5, (x,y)->linf(x,y,μ)/5);

#相空間
vec_field = scatter([0], [0], title="μ=$(μ)", labels="fixed point", color=:red, xlim=[-e,e], ylim=[-e,e], xlabel="x", ylabel="v");
for x₀ in x₀_arr
    #Runge-Kutta法
    xvec = myRK(n,x₀,f,μ);

    #プロット
    plot!(xvec[1,:], xvec[2,:], labels=false, color=:blue, linewidth=2.0);
end
quiver!(xstart_array, ystart_array,quiver=quiver_tupple, color=:black);

#保存
savefig(vec_field, "figs-2d-vectors/fig1.png");
【Juliaコード4; 負の実根が1つのとき】
#パラメータ値
μ = -2.0;

#ベクトル
xstart_array, ystart_array, quiver_tupple = make_vectors(-e, e, 0.5, (x,y)->linf(x,y,μ)/5);

#相空間
vec_field = scatter([0], [0], title="μ=$(μ)", labels="fixed point", color=:red, xlim=[-e,e], ylim=[-e,e], xlabel="x", ylabel="v");
for x₀ in x₀_arr
    #Runge-Kutta法
    xvec = myRK(n,x₀,f,μ);

    #プロット
    plot!(xvec[1,:], xvec[2,:], labels=false, color=:blue, linewidth=2.0);
end
quiver!(xstart_array, ystart_array,quiver=quiver_tupple, color=:black);

#保存
savefig(vec_field, "figs-2d-vectors/fig2.png");
【Juliaコード5; 実部が負の複素共役根のとき】
#パラメータ値
μ = -1.0;

#ベクトル
xstart_array, ystart_array, quiver_tupple = make_vectors(-e, e, 0.5, (x,y)->linf(x,y,μ)/5);

#相空間
vec_field = scatter([0], [0], title="μ=$(μ)", labels="fixed point", color=:red, xlim=[-e,e], ylim=[-e,e], xlabel="x", ylabel="v");
for x₀ in x₀_arr
    #Runge-Kutta法
    xvec = myRK(n,x₀,f,μ);

    #プロット
    plot!(xvec[1,:], xvec[2,:], labels=false, color=:blue, linewidth=2.0);
end
quiver!(xstart_array, ystart_array,quiver=quiver_tupple, color=:black);

#保存
savefig(vec_field, "figs-2d-vectors/fig3.png");
【Juliaコード6; 純虚数根のとき】
#パラメータ値
μ = 0.0;

#ベクトル
xstart_array, ystart_array, quiver_tupple = make_vectors(-e, e, 0.5, (x,y)->linf(x,y,μ)/5);

#相空間
vec_field = scatter([0], [0], title="μ=$(μ)", labels="fixed point", color=:red, xlim=[-e,e], ylim=[-e,e], xlabel="x", ylabel="v");
for x₀ in x₀_arr
    #Runge-Kutta法
    xvec = myRK(n,x₀,f,μ);

    #プロット
    plot!(xvec[1,:], xvec[2,:], labels=false, color=:blue, linewidth=2.0);
end
quiver!(xstart_array, ystart_array,quiver=quiver_tupple, color=:black);

#保存
savefig(vec_field, "figs-2d-vectors/fig4.png");
【Juliaコード7; 実部が正の複素共役根のとき】
#パラメータ値
μ = 1.0;

#ベクトル
xstart_array, ystart_array, quiver_tupple = make_vectors(-e, e, 0.5, (x,y)->linf(x,y,μ)/5);

#相空間
vec_field = scatter([0], [0], title="μ=$(μ)", labels="fixed point", color=:red, xlim=[-e,e], ylim=[-e,e], xlabel="x", ylabel="v");
for x₀ in x₀_arr
    #Runge-Kutta法
    xvec = myRK(n,x₀,f,μ);

    #プロット
    plot!(xvec[1,:], xvec[2,:], labels=false, color=:blue, linewidth=2.0);
end
quiver!(xstart_array, ystart_array,quiver=quiver_tupple, color=:black);

#保存
savefig(vec_field, "figs-2d-vectors/fig5.png");
【Juliaコード8; 正の実根が1つのとき】
#パラメータ値
μ = 2.0;

#ベクトル
xstart_array, ystart_array, quiver_tupple = make_vectors(-e, e, 0.5, (x,y)->linf(x,y,μ)/5);

#相空間
vec_field = scatter([0], [0], title="μ=$(μ)", labels="fixed point", color=:red, xlim=[-e,e], ylim=[-e,e], xlabel="x", ylabel="v");
for x₀ in x₀_arr
    #Runge-Kutta法
    xvec = myRK(n,x₀,f,μ);

    #プロット
    plot!(xvec[1,:], xvec[2,:], labels=false, color=:blue, linewidth=2.0);
end
quiver!(xstart_array, ystart_array,quiver=quiver_tupple, color=:black);

#保存
savefig(vec_field, "figs-2d-vectors/fig6.png");
【Juliaコード9; 正の実根が2つのとき】
#パラメータ値
μ = 3.0;

#ベクトル
xstart_array, ystart_array, quiver_tupple = make_vectors(-e, e, 0.5, (x,y)->linf(x,y,μ)/5);

#相空間
vec_field = scatter([0], [0], title="μ=$(μ)", labels="fixed point", color=:red, xlim=[-e,e], ylim=[-e,e], xlabel="x", ylabel="v");
for x₀ in x₀_arr
    #Runge-Kutta法
    xvec = myRK(n,x₀,f,μ);

    #プロット
    plot!(xvec[1,:], xvec[2,:], labels=false, color=:blue, linewidth=2.0);
end
quiver!(xstart_array, ystart_array,quiver=quiver_tupple, color=:black);

#保存
savefig(vec_field, "figs-2d-vectors/fig7.png");