【python】ルンゲクッタ法による数値解析(多自由度系への応用)

多自由度系でのルンゲクッタ法
- 多自由度系の運動方程式
- ルンゲクッタ法によるの運動方程式の数値解析

多自由度系でのルンゲクッタ法

前回は一自由度系の運動方程式をオイラー法とルンゲクッタ法でシミュレーションしました。
mori-memo.hateblo.jp

今回は多自由度系でも適応可能な形で実装してみたいと思います。

多自由度系の運動方程式

対象として以下のシステムの数値解析を行います。

毎度絵が雑ですいませんが、よくある外力なしの多自由度のバネーマスーダンパ系です。
質量 $m_{i}$ の台車がバネ定数 $k_{i}$ のバネと減衰係数 $c_{i}$ のダンパに $n$ 個直列につながったシステムです。この台車の変位 $x_{i}$ の時間 $t$ における運動は以下の2階の常微分方程式で表すことができます。

$i = 1$ の時
$m_{i}\dfrac{d^2 x_{i}}{dt^2} = -c_{i}\dfrac{dx_{i}}{dt} -c_{i+1}(\dfrac{dx_{i}}{dt} - \dfrac{dx_{i+1}}{dt}) - k_{i}x_{i} - k_{i+1}(x_{i} - x_{i+1})$

$2 \leq i \leq n - 1$ の時
$m_{i}\dfrac{d^2 x_{i}}{dt^2} = -c_{i}(\dfrac{dx_{i}}{dt} - \dfrac{dx_{i-1}}{dt}) -c_{i+1}(\dfrac{dx_{i}}{dt} - \dfrac{dx_{i+1}}{dt}) - k_{i}(x_{i} - x_{i-1}) - k_{i+1}(x_{i} - x_{i+1})$

$i = n$ の時
$m_{i}\dfrac{d^2 x_{i}}{dt^2} = -c_{i}(\dfrac{dx_{i}}{dt} - \dfrac{dx_{i-1}}{dt}) - k_{i}(x_{i} - x_{i-1})$

ここで、変位ベクトル $\bf{x} = \it{(x_{1}, x_{2},...,x_{n})^T}$ と速度ベクトル $\bf{v} = \it{(v_{1}, v_{2},...,v_{n})^T}$ で上記を書き直すと以下のようになります。

$\dfrac{d \bf{x}}{dt} = \bf{v}$

$\dfrac{d \bf{v}}{dt} = - \bf{M}^{-1}\bf{Cv} - \bf{M}^{-1}\bf{Kx} = \it{f}\rm{(}\bf{v}, \bf{x}\rm{)}$

$\bf{M}$ は質量マトリクス、 $\bf{C}$ は減衰マトリクス、 $\bf{K}$ は剛性マトリクスで、以下の通りです。
簡単のため、 $n=3$ の時の各マトリクスを紹介します。

$\bf{M} = \bf{I}\left( \begin{array}{c} m_{1} \\ m_{2} \\ m_{3} \end{array}\right)$
$\bf{C} = \begin{pmatrix} c_{1}+c_{2} & -c_{2} & 0 \\ -c_{2} & c_{2}+c_{3} & -c_{3} \\ 0 & c_{3} & -c_{3} \end{pmatrix}$
$\bf{K} = \begin{pmatrix} k_{1}+k_{2} & -k_{2} & 0 \\ -k_{2} & k_{2}+k_{3} & -k_{3} \\ 0 & k_{3} & -k_{3} \end{pmatrix}$

ここまでくれば、あとは前回とほぼ同様にやれば数値解析解を求めることができます。
さっそく実装してみましょう

ルンゲクッタ法によるの運動方程式の数値解析

まずは各種変数を定義します。自由度を増やしても対応可能なつくりにしたいと思います。
例として3自由度でやってみます。

# 初期値、変数の定義
m = np.array((10, 5, 10))     # 各自由度の質量 (m1, m2, m3, ... , mn)
c = np.array((1, 2, 2)) # 各自由度の減衰係数 (c1, c2, c3, ... , cn)
k = np.array((5000, 1000, 3000))    # 各自由度の弾性係数(剛性) (k1, k2, k3, ... , kn)
n = len(m)                    # 自由度
x_0 = np.array((-1.0, 0, 0))  # 変位の初期値
v_0 = np.array((0, 0, 0))     # 速度の初期値
t_max = 20                        # 終了時間
dt = 0.001                    # サンプリング間隔
t = np.arange(0, t_max, dt)       # 時間データ

ここのm, c, kを変えれば好きな自由度で解析できます。
続いて、運動方程式中の $\it{f}\rm{(}\bf{v}, \bf{x}\rm{)}$ を関数化します

# 運動方程式の関数化
def func(v, x, m, c, k):
    '''運動方程式'''
    # 各種変数定義
    n = len(m)                     # 自由度
    M = np.zeros((n, n))           # 質量マトリクス
    C = np.zeros((n, n))           # 減衰マトリクス
    K = np.zeros((n, n))           # 剛性マトリクス
    
    # マトリクスの生成
    for i in range(n):
        M[i, i] = m[i]
        if i == 0:
            C[i, [i, i + 1]] = [c[i] + c[i + 1], -c[i + 1]]
            K[i, [i, i + 1]] = [k[i] + k[i + 1], -k[i + 1]]
        elif i == n - 1:
            C[i, [-2, -1]] = [-c[i], c[i]]
            K[i, [-2, -1]] = [-k[i], k[i]]
        else:
            C[i, [i - 1, i, i + 1]] = [-c[i], c[i] + c[i + 1], -c[i + 1]]
            K[i, [i - 1, i, i + 1]] = [-k[i], k[i] + k[i + 1], -k[i + 1]]
    M_inv = np.linalg.inv(M)       # 質量マトリクスの逆行列
    return - M_inv @ C @ v.T - M_inv @ K @ x.T

@は行列の内積を行ってます。

ルンゲクッタ法も前回書いたプログラムとほぼ同じで、前回のx[i], v[i]をx[i, :], v[i, :]に書き換えてるだけです
runge_kutta_method()での返り値は変位の行列で、i行が時刻dt * iにおける各自由度の変位です。

# ルンゲクッタ法による数値解析
def runge_kutta_method(n, t, dt, x_0, v_0, func):
    '''ルンゲクッタ法'''
    x = np.zeros((len(t), n))
    v = np.zeros((len(t), n))
    for i in range(len(t) - 1):
        if i == 0:
            x[i, :] = x_0
            v[i, :] = v_0
            
        k1_x = v[i, :] * dt
        k1_v = func(v[i, :], x[i, :], m, c, k) * dt
        
        k2_x = (v[i, :] + k1_v / 2) * dt
        k2_v = func(v[i, :] + k1_v / 2, x[i, :] + k1_x / 2, m, c, k) * dt
        
        k3_x = (v[i, :] + k2_v / 2) * dt
        k3_v = func(v[i, :] + k2_v / 2, x[i, :] + k2_x / 2, m, c, k) * dt
        
        k4_x = (v[i, :] + k3_v) * dt
        k4_v = func(v[i, :] + k3_v, x[i, :] + k3_x, m, c, k) * dt
        
        x[i + 1, :] = x[i, :] + (k1_x + 2 * k2_x + 2 * k3_x + k4_x) / 6
        v[i + 1, :] = v[i, :] + (k1_v + 2 * k2_v + 2 * k3_v + k4_v) / 6
    return x

# プロット(各自由度の初期間隔を2としてプロットする)
x_rk = runge_kutta_method(n, t, dt, x_0, v_0, func)
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
for i in range(n):
    ax.plot(t, x_rk[:, i] + 2 * i, label=str(i + 1))
ax.set_xlabel('t')
ax.set_ylabel('x')
plt.legend(loc='best')
plt.show()