まえがき

25年以上前になるが、情報処理系の大学に通っていた。大学の研究室では、非線形系の解析をコンピューターを使って行っていた。当時は、FORTRANでホストコンピューターにバッチ処理を投入して、タイムシェアリング処理で計算を行ってた。あの頃は、指導されるままに行なっていたが、最近になって、内容を思い返してみて、面白い事を、やらせてもらってたのではないかと興味が湧いてきた。コンピュータ技術も発達したので、手持ちのPCでも、あの頃の研究が行えるのではないかと思い立ち、Pythonで非線形系の解析を行なってみることにした。

非線形系 Duffing-vander Pol方程式

非線形系としてDuffing-vandel Polの方程式で知られる下記の式で表される系について解析する。これは、大学時代に研究していたの時の同じ式になる。

{{d^2x}\over{dt^2}}-μ(1-γx^2){{dx}\over{dt}}+{C_1}x+x^3=Bcosνt

これを、分解して(x, y)の系に置き換える。

y={{dx}\over{dt}}

{{dy}\over{dt}}=μ(1-γx^2)y-{C_1}x-x^3+Bcosνt

更に、定数が多いので、μ=0.2, ν=1, C1=-1.0と固定した系について解析する。

非線形系の解析

上記の微分方程式をPythonを使って解く。ソースは下記になった。

import numpy as np
import matplotlib.pyplot as plt
# coding: UTF-8
from scipy.integrate import odeint, simps

def duffing(var, t, B, nu):
    """
    var = [x, p]
    dx/dt = p
    dp/dt = mu(1-gamma*x**2)*p -C1*x -x**3 + B*cos(nu*t)
    """
    mu=0.2
    gamma=1
    C1=-1.0
    
    x_dot = var[1]
    p_dot = mu*(1-gamma*var[0]**2)*var[1]-C1*var[0]-var[0]**3+B*np.cos(nu*t)

    return np.array([x_dot, p_dot])

def cfsev(var,B,nu):
    t = np.arange(0,5000,2*np.pi/nu)
    var = odeint(duffing, var, t, args=(B,nu))
    return  var.T[0][100:], var.T[1][100:]

def plot_point(var,B,nu):
    x,p=cfsev(var,B,nu)
    plt.plot(x, p, ".", markersize=4)
    plt.show()    

plot_point([1, 1],1,4)

最後の行のplot_point( [初期値], B, ν)を与えて、計算する。

一番オーソドックスなのは、安定点に落ち着くケースである。例では２点がプロットされており、νの２倍周期で安定する事を示している。

2.　この系ではリミットサイクルが観測できます。非線形系に特異な現象で、サイクル上で安定している状態である。

3.　カオス状態も存在する。カオスは、次の動作が予想できないにも関わらず、全体としては何らかの規則性を持っている特徴がある。軌道はグラフのようにカオス特有の綺麗な図形で観測できる。

このように非線形系の特徴を示す現象が観測できる興味深い系になっている。

分岐現象の観測

闇雲に変数を指定して、系の現象を確認していくのは、非効率である。そこで変数の一定範囲で変化させて現象の分布について調べてみる。加えて、系の振る舞いが変化する推移についても観測する。ソースは下記になる。

import numpy as np
import matplotlib.pyplot as plt
# coding: UTF-8
from scipy.integrate import odeint, simps

def duffing(var, t, B, nu):
    """
    var = [x, p]
    dx/dt = p
    dp/dt = mu(1-gamma*x**2)*p -C1*x -x**3 + B*cos(nu*t)
    """
    mu=0.2
    gamma=1
    C1=-1.0
    
    x_dot = var[1]
    p_dot = mu*(1-gamma*var[0]**2)*var[1]-C1*var[0]-var[0]**3+B*np.cos(nu*t)

    return np.array([x_dot, p_dot])

def cfsev(var,B,nu):
    t = np.arange(0,5000,2*np.pi/nu)
    var = odeint(duffing, var, t, args=(B,nu))
    return  var.T[0][100:], var.T[1][100:]

def orb_diagm(B,nu_min,nu_max):
    var=[0,1]
    lx=list()
    lnu=list()
    for nu in np.arange(nu_max,nu_min,-1*(nu_max-nu_min)/500):
        x,p=cfsev(var,B,nu)
        if len(x)>300:
            lx.extend(x[-300:])
            lnu.extend([nu]*300)
            var=[x[-1],p[-1]]
        print(nu);
    plt.plot(lnu, lx, ".", markersize=4)
    plt.show()
   
orb_diagm(1.0,0.3,1)

関数duffing(var, t, B, nu)とcfsev(var,B,nu)は、前回と同じである。orb_diagm(B,nu_min,nu_max)を追加している。これは、Bを固定してνをnu_minとnu_maxの間で変化させた場合の状態の推移をプロットする。最後の行のorb_diagm(1.0,0.3,1)で、Bとνを指定して実行する。

ここで、ν=2.25はカオス状態に当たる。νを2.5→2.25で見ていくと、周期4で安定していた状態から２倍に分かれて周期８になり、その後さらに分岐を繰り返してカオス状態に推移しているように見受けられる。

まとめ

Pythonを使ってDuffing-vander Pol方程式の非線形系の解析を行なってみた。安定点、リミットサイクル、カオス現象など非線形系に特有な現象が観測できた。

あとがき

大学時代は、FORTRANで大型コンピュータを使って解析していたのが、20年が経って手持ちのPCで解析できるようになった。FORTRANで何千行も書いていたコードが、Pythonのライブラリを利用する事で、20行程度で、しかも手持ちのPCで行えるようになっている。
指導を頂いていた教授からは、数年前から年賀状が来なくなった。もう80才を過ぎている御歳しになる。元気にしておられればいいが。先生、凄い世の中になりました。