grbl1.1+Arduino CNCシールドV3.5+bCNCを使用中。
BluetoothモジュールおよびbCNCのPendant機能でスマホからもワイヤレス操作可能。
その他、電子工作・プログラミング、機械学習などもやっています。
MacとUbuntuを使用。

CNCマシン全般について:
国内レーザー加工機と中国製レーザー加工機の比較
中国製レーザーダイオードについて
CNCミリングマシンとCNCルーターマシンいろいろ
その他:
利用例や付加機能など:
CNCルーター関係:



*CNCマシンの制作記録は2016/04/10〜の投稿に書いてあります。


ラベル ルンゲクッタ法 の投稿を表示しています。 すべての投稿を表示
ラベル ルンゲクッタ法 の投稿を表示しています。 すべての投稿を表示

2019年12月9日月曜日

ローレンツ・アトラクタ: オイラー法 / ルンゲクッタ法

カオス理論によく出てくるローレンツアトラクタをPython/Jupyter Notebookで描いてみました。

ローレンツ方程式:
ローレンツ方程式は、

dx/dt = -p*x + p*y
dy/dt = -x*z + r*x - y
dz/dt = x*y - b*z

であり、これら三つの微分方程式に任意の初期値x,y,zを与えると(x=0, y=0, z=0以外)、係数p=10, r=28, b=8/3の時にカオス的な振る舞いをするというもの。通常の数式なら発散/収束/周期的になりますが、そうはならないというところが面白い。また初期値がほんの少し違うだけで結果が大きく変わるという初期値鋭敏性も持つ。

初期値x=1, y=1, z=1を与えて、1ステップ0.001で50000ステップ分計算させたもの。
コードは以下。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
dt = 0.001
T = 50
p, r, b = 10, 28, 8/3
print('steps:', int(T/dt))

def lorenz_euler(V):
    x, y, z = V[0], V[1], V[2]
    X, Y, Z = [], [], []
    t = 0
    while t <= T:
        dx = (-p*x + p*y) * dt
        dy = (-x*z + r*x - y) * dt
        dz = (x*y - b*z) * dt
        x += dx
        y += dy
        z += dz
        X.append(x)
        Y.append(y)
        Z.append(z)
        t += dt
    return X, Y, Z

V0 = [1.0, 1.0, 1.0]
X, Y, Z = lorenz_euler(V0)

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection='3d')
g = 25
ax.set_xlim(-g,g)
ax.set_ylim(-g,g)
ax.set_zlim(0, g*2)
ax.plot(X, Y, Z, color='blue', alpha=0.6)
plt.show()

とりあえず、xの式に注目すると、

dx/dt = -p*x + p*y

の式を変形して、

dx = (-p*x + p*y) * dt

にすれば、dxはxにおける1ステップ分の変化量なので、whileループ内で変化量dxをxに加算していき、設定した限度を越えればwhileループを打ち切り。

x += dx

という感じで、yやzについても同様に処理し、プログラミング的には簡単な手続き。

しかし調べてみると、これは「オイラー法」という計算方法であり、あまり精度が保てないらしい。精度を出すには1ステップの刻み幅をより細かくすれば単純に向上しますが、そうすると全体のステップ数も増え計算量が多くなります。そのための工夫としてルンゲクッタ法というのがあります。以下では、各計算方法の比較をしてみます。


オイラー法:
まず、簡単な例としてオイラー法から。

f(x) = x*sin(x)

という式で試してみます。
f(x)を微分すると、

f'(x) = sin(x) + x*cos(x)

になります。
xの変化量dxは、

dx = f'(x)*dt = (sin(x) + x*cos(x))*dt

となり、

x += dx

で、xの値を更新してループ処理させます。
以下コード。
h = 0.1
T = np.pi * 2
print('steps:', int(T/h))

def f(x):
    return x * np.sin(x)

X0 = np.arange(0, T, h)
Y0 = f(X0)

def df(x): # derivative of x*sin(x) : (x*sin(x))' = sin(x) + x*cos(x) 
    return np.sin(x) + x * np.cos(x)

x = 0
y = 0
X = [x]
Y = [y]
while x <= T - h:
    dx = h
    dy = df(x) * h
    x += dx
    y += dy
    X.append(x)
    Y.append(y)

fig = plt.figure(figsize=(8, 4))
ax = fig.add_subplot(111)
ax.grid()
ax.plot(X0, Y0, label='True ', alpha=0.8, color='C0')
ax.plot(X, Y, label='Euler', alpha=0.8, color='C1')
ax.legend()
1ステップの刻み幅は0.1。合計で62ステップ。X, Yが積分したオイラー法の結果で、X0, Y0は積分を使わないで求めた結果。

わずかですが誤差があります。刻み幅をもっと細かくすれば誤差も減りますが、ローレンツアトラクタの場合は初期値鋭敏性もあるためこの誤差は致命的かもしれません。
y値最大誤差:0.3043499782


2次のルンゲクッタ法:
この方法は改良オイラー法とも呼ばれているようで、その名の通りオイラー法よりは誤差が少ない。オイラー法では刻み幅とその時のy値(高さ)を掛け合わせているため、細長い長方形を積算していますが、2次のルンゲクッタ法では長方形上端を斜め(変化量の傾き)にした台形で計算します。その分精度向上するというもの。
h = 0.1
T = np.pi * 2
print('steps:', int(T/h))

def f(x):
    return x * np.sin(x)

X0 = np.arange(0, T, h)
Y0 = f(X0)

def df(x): # derivative of sin(x)
    return np.sin(x) + x * np.cos(x)

x = 0
y = 0
X2 = [x]
Y2 = [y]

while x <= T - h:
    k1 = df(x)
    k2 = df(x + h)
    x += h
    y += h / 2 * (k1 + k2)
    X2.append(x)
    Y2.append(y)

fig = plt.figure(figsize=(8, 4))
ax = fig.add_subplot(111)
ax.grid()
ax.plot(X0, Y0, label='True ', alpha=0.8)
ax.plot(X2, Y2, label='RK2', alpha=0.8)
ax.legend()
k1とk2が台形の異なる二辺の長さ。あとは台形面積を求める計算。
今回はdtの代わりにhを刻み幅にしていますが同じ0.1です。誤差は0(h^3)程度。

かなり誤差はなくなりました。
y値最大誤差:0.0042499318


4次のルンゲクッタ法:
さらに精度をあげた方法。通常この4次のルンゲクッタ法くらいの精度が必要らしい。先ほどの2次ではk1、k2の2点でしたが、今回はk1〜k4の4点をとり、その平均を求めます。
h = 0.1
T = np.pi * 2
print('steps:', int(T/h))

def f(x):
    return x * np.sin(x)

X0 = np.arange(0, T, h)
Y0 = f(X0)

def df(x): # derivative of x*sin(x) : (x*sin(x))' = sin(x) + x*cos(x) 
    return np.sin(x) + x * np.cos(x)

x = 0
y = 0
X4 = [x]
Y4 = [y]

while x < T - h:
    k1 = df(x) * h
    k2 = df(x + h/2) * h
    k3 = df(x + h/2) * h
    k4 = df(x + h) * h
    x += h
    y += (k1 + 2*k2 + 2*k3 + k4) / 6
    X4.append(x)
    Y4.append(y)

fig = plt.figure(figsize=(8, 4))
ax = fig.add_subplot(111)
ax.grid()
ax.plot(X0, Y0, label='True ', alpha=0.8)
ax.plot(X4, Y4, label='RK4', alpha=0.8)
ax.legend()
テイラー展開により導出できるようですが、そのまま公式を覚えた方がよさそう。式中では刻み幅の両端点と中点をとって、1:2:2:1の割合で平均近似しているようです。そうすると、曲線上の中点の値df(x+h/2)と刻み幅hの両端を直線で結んだ時の中間平均値:(df(x) + df(x+h))/2には差分があるので、その2点の2/3の地点が計算で求まるようです。

グラフを見ただけでは違いはわかりにくいですが、数値的に誤差を求めてみるとさらに精度が上がっているようです。理論上、誤差は0(h^5)。
y値最大誤差:0.0000003046


scipyのライブラリ:
自前でルンゲクッタ法をコーディングするかわりに、scipyのライブラリを使ってもいいかもしれません。

scipy.integrate.solve_ivp (methodでRK45など選択可)
scipy.integrate.odeint (古いタイプ?)
scipy.integrate.RK4 (使いにくい)

solve_ivpでは計算法としてRK23やRK45も選択可能。t_evalで刻み幅を配列で指定。さらにオプションのrtol, atolで精度を指定することもでできるようなので、アルゴリズム的な工夫はほぼ必要なし。デフォルトではrtol=1e-3、atol=1e-6になっており、そのままでも大丈夫そう。


RK2、RK4、Odeint、Solve_ivpの比較:
ルンゲクッタ法(RK2、RK4)、ならびにライブラリのOdeint、Solver_ivp(RK45)をローレンツアトラクで比較してみることに(とは言っても、初期値鋭敏性のあるローレンツアトラクタなので、精度の比較はしにくいですが)。
さらに、3D/2D表示、XYZの各値の変化を1D表示するファンクションも追加しておきました。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.integrate import odeint
from scipy.integrate import solve_ivp

%matplotlib inline
#%matplotlib notebook

def df(t, xyz):
    x, y, z = xyz
    return np.array([-p*x + p*y, -x*z + r*x - y, x*y - b*z])

def lorenz_RK2(V):
    xyz = V
    t = 0
    X, Y, Z = [], [], []
    while t <= T:
        k1 = df(t, xyz)
        k2 = df(t+h, xyz + h*k1)
        xyz += (k1 + k2) * h / 2
        X.append(xyz[0])
        Y.append(xyz[1])
        Z.append(xyz[2])
        t += h
    return X, Y, Z

def lorenz_RK4(V):
    xyz = V
    t = 0
    X, Y, Z = [], [], []
    while t <= T:
        k1 = df(t, xyz)
        k2 = df(t+h/2, xyz + h*k1/2)
        k3 = df(t+h/2, xyz + h*k2/2)
        k4 = df(t+h, xyz + h*k3)
        xyz += (k1 + 2*k2 + 2*k3 + k4) * h / 6
        X.append(xyz[0])
        Y.append(xyz[1])
        Z.append(xyz[2])
        t += h
    return X, Y, Z

def lorenz_Odeint(V):
    xyz = V
    t = np.arange(0, T, h)
    X, Y, Z = odeint(df, xyz, t, tfirst=True).T 
    #X, Y, Z = odeint(f2, xyz, t, rtol=1e-9, atol=1e-9).T
    return X, Y, Z

def lorenz_ivp(V):
    xyz_0 = np.array(V)
    t = np.arange(0, T, h) 
    res = solve_ivp(df, (0, T), xyz_0, method='RK45', t_eval=t)#, rtol=1e-7, atol=1e-9)
    X = res.y[0]
    Y = res.y[1]
    Z = res.y[2]
    print(X[-1], Y[-1], Z[-1])
    return X, Y, Z

def plot3D(methods, labels):
    from matplotlib.colors import BASE_COLORS
    colorlist = list(BASE_COLORS.keys())
    
    fig = plt.figure(figsize=(10,10))
    ax = fig.add_subplot(111, projection='3d')
    g = 30 
    ax.set_xlim(-g, g)
    ax.set_ylim(-g, g)
    ax.set_zlim(0, 2*g)
    ax.xaxis.pane.set_facecolor('white')
    ax.yaxis.pane.set_facecolor('white')
    ax.zaxis.pane.set_facecolor('white')

    for i, XYZ in enumerate(methods):
        X = XYZ[0]
        Y = XYZ[1]
        Z = XYZ[2]
        
        #color = 'C' + str(i)
        color = colorlist[i]
        ax.plot(X, Y, Z, color=color, alpha=0.5, label=labels[i]) # plot method
        ax.scatter(X[0], Y[0], Z[0], color=color, s=60, marker='x') # start point
        ax.text(X[0], Y[0], Z[0], s='START  ', color=color, verticalalignment='top', horizontalalignment='right')
        ax.scatter(X[-1], Y[-1], Z[-1], color=color, s=30, marker='o') # end point

    plt.legend()
    plt.show()

def plot2D(methods, labels, axis):
    from matplotlib.colors import BASE_COLORS
    colorlist = list(BASE_COLORS.keys())
    
    fig = plt.figure(figsize=(8,8))
    ax = fig.add_subplot(111)
    ax.axis([-30, 30, 0, 60], 'equal')
    if axis == 'XY':
        ax.axis([-30, 30, -30, 30], 'equal')     
    ax.grid()
    
    for i, XYZ in enumerate(methods):
        if axis == 'XY':
            A = XYZ[0]
            B = XYZ[1]
        elif axis == 'XZ':
            A = XYZ[0]
            B = XYZ[2]
        else:
            A = XYZ[1]
            B = XYZ[2]
        
        color = colorlist[i]
        ax.plot(A, B, color=color, alpha=0.5, label=labels[i])  # plot line
        ax.scatter(A[0], B[0], color=color, s=60, marker='x')   # start point
        ax.text(A[0], B[0], s='START  ', color=color, verticalalignment='center', horizontalalignment='right')
        ax.scatter(A[-1], B[-1], color=color, s=30, marker='o') # end point
    plt.legend()
    plt.show()
    print(axis + ' coordinates')

def plot1D(methods, labels, axis):
    plt.figure(figsize=(14, 3))
    for i, XYZ in enumerate(methods):
        if axis == 'X':
            A = XYZ[0]
        elif axis == 'Y':
            A = XYZ[1]
        else:
            A = XYZ[2]
        plt.plot(list(range(len(A))), A, color='C'+str(i), alpha=0.6, label=labels[i])
    plt.legend()
    plt.show()

実行用のコードは以下。
h = 1e-3 # dt
T = 30   # total time
p, r, b = 10, 28, 8/3
print('steps:', int(T/h))

xyz_0 = [10.0, 10.0, 10.0] # initial values for the methods 
methods = [lorenz_RK2(xyz_0), lorenz_RK4(xyz_0), lorenz_Odeint(xyz_0), lorenz_ivp(xyz_0)]
labels = ['RK2', 'RK4', 'ODE', 'IVP']

plot3D(methods, labels)
# plot2D(methods, labels, 'XY')
# plot2D(methods, labels, 'XZ')
plot2D(methods, labels, 'YZ')
plot1D(methods, labels, 'X')

表示結果。
 4種類の異なる計算方法。x印が開始点、●印が終着点(水色のIVPは見えにくくなっていますが)。


YZ座標。この方向からみるとほぼ線対称の形になります。


縦軸はX値、横軸はステップ数。
6000ステップくらいまでは一緒ですが、その後は次第にずれてきています。単体だけ見ると、規則的に見える部分もありますが、やはり非周期的になっています。


初期値鋭敏性:
計算方法は同じで、初期値をほんの少しだけ変えて実験してみました。
初期値A:X=10.0、Y=10.0、Z=10.0
初期値B:X=10.001、Y=10.0、Z=10.0
BのXの初期値だけ0.001違います。あとは同じ条件で50000ステップ。
以下が実行コード。
h = 1e-3
T = 50
p, r, b = 10, 28, 8/3
print('steps:', int(T/h))

xyz_0 = [10.0, 10.0, 10.0]
xyz_1 = [10.001, 10.0, 10.0]
methods = [lorenz_RK4(xyz_0), lorenz_RK4(xyz_1)]
labels = ['RK4:'+str(xyz_0), 'RK4:'+str(xyz_1)]

plot3D(methods, labels)
# plot2D(methods, labels, 'XY')
# plot2D(methods, labels, 'XZ')
# plot2D(methods, labels, 'YZ')
plot1D(methods, labels, 'X')
plot1D(methods, labels, 'Z')

結果表示。
軌道に差があるため、青と緑のムラができています。✖️が開始点、●は終着点。


上が50000ステップにおけるX値。下はZ値。
10000ステップを超えたあたりから差がでてきます。ずれるというより、異なったパターンになってきています。
完璧なランダムとは言えませんが、数多くのステップ数を繰り返しても周期的にならないという部分が興味深い。


分布について:
ランダムの指標として、値全体の分布もとってみました。
これはオイラー法による結果ですが、上から、X値、Y値、Z値の分布です(50000ステップ)。
基本的には中央に寄っていますが、正規分布とも違う。


まとめ:
元々はカオス的な振る舞いをするローレンツアトラクタを利用して乱数生成できないかということを考えておりましたが、計算方法であるルンゲクッタ法やライブラリについての実験となりました。これはこれでいろいろ勉強になります。精度を高める計算方法にはいくつかあるようですが、PID制御もそのうちの一つという感じです。ということで、引き続き調査してみたいと思います。




人気の投稿