ノートブック

古き良きスタイルのノートブックです。

シンクロするモデルのメモ

リミットサイクルのモデルと簡単な同期するモデルをいくつかメモしておきます。

ファン・デル・ポール オシレーター

一般的なかたちは以下のように表され

\begin{align*} \dot{x} = v \tag{1-1} \end{align*}

\begin{align*} \dot{v} = \mu (1-x^{2})v - x \tag{1-2} \end{align*}

パラメータを適当に設定して、いくつかの初期条件からの挙動を以下のスクリプトでプロットしてみます。

# -*- coding: utf-8 -*-

import numpy as np
import time,sys,math,random
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

def van_der_pol():

    F = 0
    x = random.uniform(-3.0,3.0)
    v = random.uniform(-3.0,3.0)
    dt = 0.01
    mu = 1.0
    plot_x = []
    plot_v = []
    
    for i in range(1000):
        F = mu*(1.0-x*x)*v-x
        v = v + F*dt
        x = x + v*dt
        plot_x.append(x)
        plot_v.append(v)

    plt.axes().set_aspect('equal', 'datalim')
    plt.xlabel("x")
    plt.ylabel("v")
    plt.plot(plot_x, plot_v, color="black", linewidth=0.5, linestyle="-")
    
    #time series plotting
    #t = np.arange(0.0, 100.0, 0.01)
    #plt.plot(t, plot_x, color="black", linewidth=0.5, linestyle="-")
    #plt.axis([0, 100, -10, 10])

for i in range(20):
    van_der_pol()
plt.show()

プロットは以下のように出力されます。

f:id:primitivem:20171119224214p:plain

別のリミットサイクル オシレーター

ロボットを制御するCPGの要素となるリミットサイクルを探していてこのレビュー論文っぽいものを見つけた

Auke Jan Ijspeert (2008). Central pattern generators for locomotion control in animals and robots: A review. Neural Networks, 21, 642–653

ので、ここに載っている円形の位相空間を示すモデルを自分の手元で動かしてみました。
方程式は以下です。

 {
\displaystyle
\begin{align}
\tau \dot{x} = v \tag{2-1}
\end{align}
}

 {
\displaystyle
\begin{align}
\tau \dot{v} = - \alpha \frac{x^{2} + v^{2} - E}{E} v - x \tag{2-2}
\end{align}
}

これを適当なパラメータ設定の下で以下のスクリプトで解を描いてみます。

import numpy as np
import time,sys,math,random
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

def limit_cycle():

    F = 0
    x = random.uniform(-1.5,1.5)
    v = random.uniform(-1.5,1.5)
    dt = 0.01
    tau = 1.0
    alpha = 1.0
    E = 1.0
    plot_x = []
    plot_v = []
    
    for i in range(10000):
        F = ( -alpha*(x*x+v*v-E)*v/E - x )/tau
        v = v + F*dt
        x = x + v*dt
        
        plot_x.append(x)
        plot_v.append(v)

    plt.axes().set_aspect('equal', 'datalim')
    plt.xlabel("x")
    plt.ylabel("v")
    plt.plot(plot_x, plot_v, color="black", linewidth=0.5, linestyle="-")

    #time series plotting
    #t = np.arange(0.0, 100.0, 0.01)
    #plt.plot(t, plot_x, color="black", linewidth=0.5, linestyle="-")
    #plt.axis([0, 100, -10, 10])

for i in range(20):
    limit_cycle()
plt.show()

出力はこんな感じ。

f:id:primitivem:20171119225138p:plain

クラモト モデル

単純な周期運動をする振動子を接続することで同期現象が発生する、いわゆるクラモト・モデルも動かして試してみます。 式は以下のようなもので

\begin{align*} \frac{d\theta_{i}}{dt} = \omega_{i} + \frac{K}{N} \sum_{j=1}^{N} sin(\theta_{j}-\theta_{i}) \tag{3-1} \end{align*}

これの挙動を以下のスクリプトで確かめてみます。

import numpy as np
import time,sys,math,random
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

def kuramoto_model_gif():

    N = 25
    dt = 0.01
    theta = []
    omega = []
    K = 1.0
    for i in range(N):
        theta.append(random.uniform(0,6.28))
        omega.append(10.0)

    #prepare plotting field
    fig = plt.figure(figsize=(6.2,6))
    plt.axes(xlim=(0, 100), ylim=(0, 100))
    circle = []
    for i in range(N):
        x = 10+20*(i%5)
        y = 10+20*(i/5)
        circle.append( plt.Circle( (x, y), radius=1, fc='black') )
        plt.gca().add_patch(circle[i])

    def update(frame_number):
        for step in range(5):
            for i in range(N):
                connect = 0.0
                for j in range(N):
                    connect += math.sin(theta[i]-theta[j])

                theta[i] = theta[i] + (omega[i] - (K/float(N))*connect)*dt
                if theta[i] > 6.28:
                    theta[i] = 0.0

        for i in range(N):
            circle[i].radius = theta[i]
            
    animation = FuncAnimation(fig, update, interval=1, frames=120)
    animation.save("output.gif", writer="imagemagick")
    #plt.show()

kuramoto_model_gif()

アニメーションをGIFとして出力されたものが以下です。

f:id:primitivem:20171119230844g:plain

以上簡単なシンクロするモデルのメモでした。