微分方程式の数値計算 線形加速度法(ニューマークβ法)

線形加速度法概要

　振動の問題によく用いられる解法に線形加速度法というものがある。振動の変位$u(t+Δt)$と速度$\dot{u}(t+Δt)$以下の式で計算する。

\begin{align}
u(t+Δt) &= u(t)+ Δt\dot{u}(t)+\frac{(Δt)^2}{2!}\ddot{u}(t)+\frac{(Δt)^3}{3!}\frac{\ddot{u}(t+Δt)-\ddot{u}(t)}{Δt}\\
\dot{u}(t+Δt) &= \dot{u}(t)+Δt\frac{\ddot{u}(t)+\ddot{u}(t+Δt)}{2}\\
\end{align}

変位$u(t+Δt)$はテイラー展開の3次までの項に似ているが、3次の項は加速度が１次式で変化する(線形加速度)と仮定した式になっている。速度$\dot{u}(t+Δt)$についても同様で、いわゆる台形公式である。

陽解法と陰解法

　オイラー法やルンゲクッタ法のような数値的計算方法は「陽解法」と呼ばれる。一方、時刻$(t+Δt)$時点の計算をするために、時刻$(t+Δt)$時点の値を使って計算するような数値的計算方法は「陰解法」と呼ばれる。線形加速度法は陰解法である。
陰解法の計算では一般的には方程式を満たす値を収束演算にて逐次微修正しながら解を求めていくことになるが、計算安定性は陽解法に比べて安定している。

計算安定性と計算精度について

　計算安定性と計算精度は異なることに注意する。ルンゲクッタ法は4次精度であるが線形加速度法はそれ未満である。両数値計算手法で解けるような計算ステップサイズであればルンゲクッタ法の計算結果の方が真値に近い値かもしれないが、計算安定性が満たされない計算ステップサイズでルンゲクッタ法の計算が発散してしまい、線形加速度法の計算結果は発散せずそれなりの計算結果を示した場合は、線形加速度法の精度が良いのではなく、あくまで安定性が良いということである。発散はしないが真値とかけ離れていることはあり得る。シミュレーション結果と各種数値計算法の特徴から結果の妥当性を吟味することは解析者にとっては必要である。

MCK型(2階微分型)微分方程式の数値計算の狙い

　陰解法は収束演算が伴うため計算コストが高い。しかし、陽解法のように代数的に解くことができれば陰解法の計算安定性が良いという利点と計算の簡易さの両立が可能となる。そのようなテクニックがないかを考える。
今、未知変数が$u(t+Δt)、\dot{u}(t+Δt)、\ddot{u}(t+Δt)$に対して式が2個しかないため、もう一つ追加することを考える。そこで、今回対象とする物理現象の運動方程式を導入する。質量$m$、減衰$c$、剛性$k$として、以下のように時刻$(t+Δt)$時点で満たすべき運動方程式を立てる。

m\ddot{u}(t+Δt)+c\dot{u}(t+Δt)+ku(t+Δt) = f(t+Δt)

この式を用いて3式の連立方程式より$u(t+Δt)、\dot{u}(t+Δt)、\ddot{u}(t+Δt)$を求めることで、陰解法を代数的に解く。($f(t+Δt)$は入力条件として既知であるとする。)
線形加速度法(ニューマークβ法)、ウィルソンθ法やフーボルト法は陰解法を代数的に解く数値計算手法で、上記のようにmckの2階微分方程式を利用することからMCK型(2階微分型)と呼ばれることがある。

※以下の流体分野の参考書籍では流体の性質を利用した例「半陰解法」が掲載されている。
越塚誠一著「粒子法入門 - 流体シミュレーションの基礎から並列計算と可視化まで」丸善出版 2014.

線形加速度法の直接計算法

以下の流れで解を求める。
①$u(t+Δt)、\dot{u}(t+Δt)$を運動方程式に代入する。
②$\ddot{u}(t+Δt)$について解く。
③得られた$\ddot{u}(t+Δt)$を$u(t+Δt)、\dot{u}(t+Δt)$の式に代入し解を求める。

具体的計算(1自由度の場合)

$u(t+Δt)、\dot{u}(t+Δt)$を運動方程式に代入すると以下のようになる。

\begin{align}
&m\ddot{u}(t+Δt)+c\left(\dot{u}(t)+Δt\frac{\ddot{u}(t)+\ddot{u}(t+Δt)}{2}\right) \\
&+k\left(u(t)+ Δt\dot{u}(t)+\frac{(Δt)^2}{2!}\ddot{u}(t)+\frac{(Δt)^3}{6}\frac{\ddot{u}(t+Δt)-\ddot{u}(t)}{Δt}\right) = f(t+Δt) \\
\end{align}

少し整理すると以下のようになる。

\begin{align}
&m\ddot{u}(t+Δt)+c\left(\dot{u}(t)+Δt\frac{\ddot{u}(t)+\ddot{u}(t+Δt)}{2}\right) \\
&+k\left(u(t)+ Δt\dot{u}(t)+\frac{(Δt)^2}{3}\ddot{u}(t)+\frac{(Δt)^2}{6}\ddot{u}(t+Δt)\right) = f(t+Δt)
\end{align}

$\ddot{u}(t+Δt)$について解くと以下のようになる。

\therefore \ddot{u}(t+Δt)=\frac{f(t+Δt)-c\left(\dot{u}(t)+\frac{Δt}{2}\ddot{u}(t)\right) - k\left(u(t)+ Δt\dot{u}(t)+\frac{(Δt)^2}{3}\ddot{u}(t)\right)}{m+\frac{Δt}{2}c+\frac{(Δt)^2}{6}k}

得られた$\ddot{u}(t+Δt)$を、以下の$u(t+Δt)、\dot{u}(t+Δt)$の式に代入し解を求める。

\begin{align}
\dot{u}(t+Δt) &= \dot{u}(t)+Δt\frac{\ddot{u}(t)+\ddot{u}(t+Δt)}{2}\\
u(t+Δt) &= u(t)+ Δt\dot{u}(t)+\frac{(Δt)^2}{3}\ddot{u}(t)+\frac{(Δt)^2}{6}\ddot{u}(t+Δt)\\
\end{align}

※最初に求める変数は、$\dot{u}(t+Δt)$の場合、$u(t+Δt)$の場合でも良い。

具体的計算(多自由度の場合)

$u(t+Δt)、\dot{u}(t+Δt)$、$\ddot{u}(t+Δt)$が多自由度の場合でも計算可能である。その場合は行列演算になることに注意して以下のようになる。(質量マトリクス$M$、減衰マトリクス$C$、剛性マトリクス$K$とする。)

\begin{align}
\vec{\ddot{u}}(t+Δt)&= \left(M+\frac{Δt}{2}C+\frac{(Δt)^2}{6}K\right)^{-1} \left(\vec{f}(t+Δt)-C\left(\vec{\dot{u}}(t)+\frac{Δt}{2}\vec{\ddot{u}}(t)\right) - K\left(\vec{u}(t)+ Δt\vec{\dot{u}}(t)+\frac{(Δt)^2}{3}\vec{\ddot{u}}(t)\right)\right)\\
\vec{\dot{u}}(t+Δt) &= \vec{\dot{u}}(t)+Δt\frac{\vec{\ddot{u}}(t)+\vec{\ddot{u}}(t+Δt)}{2}\\
\vec{u}(t+Δt) &= \vec{u}(t)+ Δt\vec{\dot{u}}(t)+\frac{(Δt)^2}{3}\vec{\ddot{u}}(t)+\frac{(Δt)^2}{6}\vec{\ddot{u}}(t+Δt)\\
\end{align}

ニューマークβ法

　線形加速度法を一般化した公式として、ニューマーク$\beta$法がある。これは、線形加速度法の3次項の係数をβとした以下のような式である。

\begin{align}
u(t+Δt) &= u(t)+ Δt\dot{u}(t)+\frac{(Δt)^2}{2!}\ddot{u}(t)+\beta(Δt)^3\frac{\ddot{u}(t+Δt)-\ddot{u}(t)}{Δt}\\
&= u(t)+ Δt\dot{u}(t)+\frac{(Δt)^2}{2!}\ddot{u}(t)+\beta(Δt)^2 \left( \ddot{u}(t+Δt)-\ddot{u}(t) \right)\\
\dot{u}(t+Δt) &= \dot{u}(t)+Δt\frac{\ddot{u}(t)+\ddot{u}(t+Δt)}{2}\\
\end{align}

$\beta$ は以下の範囲とする。

0 \leq \beta \leq \frac{1}{2}

$\beta$ の設定により公式の特徴が変わる。

β=0 の場合：

　u(t+Δt) は以下のようになり、陽解法に近い性質となる。

u(t+Δt) = u(t)+ Δt\dot{u}(t)+\frac{(Δt)^2}{2!}\ddot{u}(t)

β=1/2 の場合：

　u(t+Δt) は以下のようになり、陰解法に近い性質となる。

u(t+Δt) = u(t)+ Δt\dot{u}(t)+\frac{(Δt)^2}{2!}\ddot{u}(t+Δt)

β=1/6 の場合：

　線形加速度法の式になる。多くの場合、β=1/6としているため、ニューマークβ法は線形加速度法の別名のようなものである。

β=1/４ の場合：

　u(t+Δt) は以下のようになる。

\begin{align}
u(t+Δt) &= u(t)+ Δt\dot{u}(t)+\frac{(Δt)^2}{2!}\ddot{u}(t)+\frac{1}{4}(Δt)^2 \left( \ddot{u}(t+Δt)-\ddot{u}(t) \right)\\
 &= u(t)+ Δt\dot{u}(t)+ \left( \frac{(Δt)^2}{2}\ddot{u}(t)-\frac{(Δt)^2 }{4}\ddot{u}(t) \right) + \frac{(Δt)^2 }{4}\ddot{u}(t+Δt)\\
&= u(t)+ Δt\dot{u}(t)+ \frac{(Δt)^2}{4}\ddot{u}(t) + \frac{(Δt)^2 }{4}\ddot{u}(t+Δt)\\
&= u(t)+ Δt\dot{u}(t)+ \frac{(Δt)^2}{2!} \frac{\ddot{u}(t) + \ddot{u}(t+Δt)}{2}\\
\end{align}

上記の式は、加速度を時刻tから(t+Δt)までの平均$(\ddot{u}(t) + \ddot{u}(t+Δt))/2$ で一定値加速しているとみなした式となっており、「平均加速度法」と呼ばれている。

ニューマークβ法の直接計算法

線形加速度法と同様の流れで解を求める。
①$u(t+Δt)、\dot{u}(t+Δt)$を運動方程式に代入する。
②$\ddot{u}(t+Δt)$について解く。
③得られた$\ddot{u}(t+Δt)$を$u(t+Δt)、\dot{u}(t+Δt)$の式に代入し解を求める。

具体的計算(1自由度の場合)

$u(t+Δt)、\dot{u}(t+Δt)$を運動方程式に代入すると以下のようになる。

\begin{align}
&m\ddot{u}(t+Δt)+c\left(\dot{u}(t)+Δt\frac{\ddot{u}(t)+\ddot{u}(t+Δt)}{2}\right) \\
&+k\left( u(t)+ Δt\dot{u}(t)+\frac{(Δt)^2}{2!}\ddot{u}(t)+\beta(Δt)^2 \left( \ddot{u}(t+Δt)-\ddot{u}(t) \right) \right) = f(t+Δt) \\
\end{align}

$\ddot{u}(t+Δt)$について解くと以下のようになる。

\therefore \ddot{u}(t+Δt)=\frac{f(t+Δt)-c\left(\dot{u}(t)+\frac{Δt}{2}\ddot{u}(t)\right) - k\left(u(t)+ Δt\dot{u}(t)+(\frac{1}{2}-\beta)(Δt)^2\ddot{u}(t)\right)}{m+\frac{Δt}{2}c+\beta(Δt)^2k}

得られた$\ddot{u}(t+Δt)$を、以下の$u(t+Δt)、\dot{u}(t+Δt)$の式に代入し解を求める。

\begin{align}
\dot{u}(t+Δt) &= \dot{u}(t)+Δt\frac{\ddot{u}(t)+\ddot{u}(t+Δt)}{2}\\
u(t+Δt) &= u(t)+ Δt\dot{u}(t)+\frac{(Δt)^2}{2}\ddot{u}(t)+\beta(Δt)^2 \left( \ddot{u}(t+Δt)-\ddot{u}(t) \right)
\end{align}

具体的計算(多自由度の場合)

※こちらも線形加速度法同様のため割愛

Pythonによるシミュレーション確認

　Pythonによる線形加速度法(ニューマークβ法)の動作を計算シミュレーションで確認してみる。なお、解析解と比較できるよう簡単な1次元の振動問題を取り上げる。
(こちらのオイラー法やルンゲクッタ法で取り上げた例題と同様の問題設定である。)

問題設定

以下のような1自由度の振動モデルを考える。$m$を質点の質量、$c$を減衰、$k$をバネ定数とし、$x$が質点座標を表す。$\dot{x}$は質点速度、$\ddot{x}$は質点の加速度を表し、初期変位$x_0$と初速度$v_0$を与える。

上式の運動方程式は以下となる。

m\ddot{x}+c\dot{x}+kx = 0

解析解

解析解は以下のようになる。(導出過程はこちらオイラー法やルンゲクッタ法を参照)

\begin{align}
x &= Ae^{-\zeta\omega _nt}cos(\omega _dt - \phi) \\
A &= \sqrt{x_0^2 + \Bigl(\frac{\zeta\omega _nx_0+v_0}{\omega _d} \Bigr)^2} \\
\phi &= tan^{-1}\Bigl(\frac{\zeta\omega _nx_0+v_0}{x_0 \omega _d}\Bigr)
\end{align}

作成したPythonプログラム

import numpy as np
import matplotlib.pyplot as plt
import sympy as sp

cSimSettings = {
    "SIM_START": 0.0,
    "SIM_END": 5.0,
    "SIM_STEP": 0.01,
}

cMdlParams = {
    "M": 5.0,
    "C": 16.0,
    "K": 320.0,
    "X0": 0.05,         # 初期変位
    "V0": 0.4,          # 初期速度
    "A0": 0.0,
    "F": 0.0,
}

cNewmarkParams = {
    "BETA": 1.0/6.0
}


def calcAnalyticalSolution(t):
    # モデルパラメータ設定
    m = cMdlParams["M"]
    c = cMdlParams["C"]
    k = cMdlParams["K"]
    v0 = cMdlParams["V0"]
    x0 = cMdlParams["X0"]
    # 振動パラメータ計算
    wn = np.sqrt(k/m)                            # 不減衰固有角振動数
    zeta = c/(2.0*np.sqrt(m*k))                  # 減衰比
    wd = np.sqrt(1-zeta**2)*wn                   # 減衰固有角振動数
    # 解の計算
    Amp = np.sqrt((x0**2 + ((zeta*wn*x0 + v0)/wd)**2))
    phi = np.arctan((zeta*wn*x0 + v0)/(x0*wd))
    x = Amp*np.exp(-zeta*wn*t)*np.cos(wd*t - phi)
    return x


def definiteSolverNewmarkBeta():
    # 代数計算用シンボル定義
    m, c, k = sp.symbols('m, c, k')
    an, vn, xn = sp.symbols('an, vn, xn')
    an1, vn1, xn1, fn1 = sp.symbols('an1, vn1, xn1, fn1')
    dt, beta = sp.symbols('dt, beta')

    # (t+Δt)時点で満たすべき1次元振動モデルの運動方程式の定義
    eqm = m*an1 + c*vn1 + k*xn1 - fn1
    eqm = sp.lambdify((vn1, xn1), eqm, "numpy")  # 速度と位置を変数化

    # 速度と位置の数値積分式の定義
    vn1 = vn + dt*(an + an1)/2
    xn1 = xn + dt*vn + dt**2/2*an + beta*dt**2*(an1 - an)

    # 運動方程式に代入して加速度について解く
    fan1 = eqm(vn1, xn1)
    fan1 = sp.solve(fan1, an1)[0]  # 配列形式となって返ってきたため要素を取り出すために[0]を追加

    # 加速度、速度、位置の計算関数化
    an1_solver = sp.lambdify((m,c,k,xn,vn,an,fn1,dt,beta), fan1, "numpy")
    vn1_solver = sp.lambdify((vn,an,an1,dt), vn1, "numpy")
    xn1_solver = sp.lambdify((xn,vn,an,an1,dt,beta), xn1, "numpy")
    
    return(an1_solver, vn1_solver, xn1_solver)


def solverNewmarkBeta(m,c,k,xn,vn,an,fn1,dt,beta):
    an1_num = fn1 - c*(vn + dt*an/2.0) - k*(xn + dt*vn + (1.0/2.0 - beta)*(dt**2)*an)
    an1_den = m + dt*c/2.0 + beta*(dt**2)*k
    an1 = an1_num/an1_den
    vn1 = vn + dt*(an + an1)/2.0
    xn1 = xn + dt*vn + (dt**2)*an/2.0 + beta*(dt**2)*(an1 - an)
    return(an1, vn1, xn1)


if __name__ == "__main__":
    # 変数準備(固定ステップで計算するためあらかじめ用意しておく)
    sim_start = cSimSettings["SIM_START"]
    sim_end = cSimSettings["SIM_END"]
    sim_step = cSimSettings["SIM_STEP"]
    num_step = int((sim_end-sim_start)/sim_step) + 1
    t = np.linspace(sim_start, sim_end, num_step)

    # 数値計算の格納結果
    a_nms = np.zeros(num_step)
    v_nms = np.zeros(num_step)
    x_nms = np.zeros(num_step)
    a_nm = np.zeros(num_step)
    v_nm = np.zeros(num_step)
    x_nm = np.zeros(num_step)

    # 解析解の計算
    num_step = int((sim_end-sim_start)/0.01) + 1
    ta = np.linspace(sim_start, sim_end, num_step)
    x = calcAnalyticalSolution(ta)

    # 数値計算による解
    a0 = cMdlParams["A0"]
    v0 = cMdlParams["V0"]
    x0 = cMdlParams["X0"]
    m = cMdlParams["M"]
    c = cMdlParams["C"]
    k = cMdlParams["K"]
    f = cMdlParams["F"]
    dt = sim_step
    b = cNewmarkParams["BETA"]
    calc_an1, calc_vn1, calc_xn1 = definiteSolverNewmarkBeta()
    for i in range(t.size):
        if (i < 1):
            # 初期値
            [a_nms[i], v_nms[i], x_nms[i]] = [a0, v0, x0]
            [a_nm[i], v_nm[i], x_nm[i]] = [a0, v0, x0]
        else:
            # シンボル計算で作成したニューマークβ法による計算
            a_nms[i] = calc_an1(m,c,k,x_nms[i-1],v_nms[i-1],a_nms[i-1],f,dt,b)
            v_nms[i] = calc_vn1(v_nms[i-1],a_nms[i-1],a_nms[i],dt)
            x_nms[i] = calc_xn1(x_nms[i-1],v_nms[i-1],a_nms[i-1],a_nms[i],dt,b)
            # ニューマークβ法による計算
            a_nm[i],v_nm[i], x_nm[i] = solverNewmarkBeta(m,c,k,
                                                         x_nm[i-1],
                                                         v_nm[i-1],
                                                         a_nm[i-1],
                                                         f,dt,b)

    # グラフ確認
    plt.plot(ta, x, 'k', lw=3.0, label='Analytical Solution')
    plt.plot(t, x_nms, 'g', lw=2.0, label='Newmark-β method(Symbolic)')
    plt.plot(t, x_nm, 'r--', lw=1.0, label='Newmark-β method')
    plt.xlim(0, 3.5)
    plt.ylim(-0.04, 0.07)
    plt.xlabel('Time (s)')
    plt.ylabel('Displacement x(m)')
    plt.legend()
    plt.grid()
    plt.show()
    #plt.savefig("./1dofVibSim_Newmarkbeta.png")

解析解は関数「calcAnalyticalSolution」で計算している。
線形加速度法(ニューマークβ法)は今回2通りの実装方法（「definiteSolverNewmarkBeta」と「solverNewmarkBeta」）を試した。

関数 definiteSolverNewmarkBeta とそれを用いた数値計算の簡単な説明

こちらはscipyのシンボリック計算を利用して線形加速度法(ニューマークβ法)の計算式を導出している。
関数冒頭でまず計算に使用するシンボルを定義している。

    # (t+Δt)時点で満たすべき1次元振動モデルの運動方程式の定義
    eqm = m*an1 + c*vn1 + k*xn1 - fn1
    eqm = sp.lambdify((vn1, xn1), eqm, "numpy")  # 速度と位置を変数化

こちらで、時刻t+Δt時点で満たす運動方程式を定義している。各種変数はそれぞれ以下を表している。
an1:$\ddot{u}(t+Δ)$
vn1:$\dot{u}(t+Δ)$
xn1:$u(t+Δ)$
fn1:$f(t+Δ)$
上記2行目にて、線形加速度法(ニューマークβ法)の直接計算法にて解を求める流れ①$u(t+Δ)$、$\dot{u}(t+Δ)$を運動方程式に代入できるように、引数で受け取るufunc関数eqm を定義している。

    # 速度と位置の数値積分式の定義
    vn1 = vn + dt*(an + an1)/2
    xn1 = xn + dt*vn + dt**2/2*an + beta*dt**2*(an1 - an)

    # 運動方程式に代入して加速度について解く
    fan1 = eqm(vn1, xn1)
    fan1 = sp.solve(fan1, an1)[0]  # 配列形式となって返ってきたため要素を取り出すために[0]を追加

こちらで$u(t+Δ)$、$\dot{u}(t+Δ)$を定義し運動方程式に代入している。上記の最終行にて、解を求める流れ②$\ddot{u}(t+Δ)$について解いている。sp.solve(fan1, an1) で実行した際、fan1に配列形式となって返ってくることにより後の計算でエラーが出たため、要素を取り出すために[0]を追加している。

    # 加速度、速度、位置の計算関数化
    an1_solver = sp.lambdify((m,c,k,xn,vn,an,fn1,dt,beta), fan1, "numpy")
    vn1_solver = sp.lambdify((vn,an,an1,dt), vn1, "numpy")
    xn1_solver = sp.lambdify((xn,vn,an,an1,dt,beta), xn1, "numpy")

各変数の計算式が出来上がったため関数化している。
main関数ではまずこれらの関数が使えるように以下を実行する。

calc_an1, calc_vn1, calc_xn1 = definiteSolverNewmarkBeta()

各時刻での計算は以下のように行っている。

# シンボル計算で作成したニューマークβ法による計算
a_nms[i] = calc_an1(m,c,k,x_nms[i-1],v_nms[i-1],a_nms[i-1],f,dt,b)
v_nms[i] = calc_vn1(v_nms[i-1],a_nms[i-1],a_nms[i],dt)
x_nms[i] = calc_xn1(x_nms[i-1],v_nms[i-1],a_nms[i-1],a_nms[i],dt,b)

プログラム内での計算じっくステップi は、説明中の時刻$t+Δt$に相当する。まず加速度について解き、解を求める流れ③のように求めた加速度a_nms[i]を使って速度v_nms[i]と位置x_nms[i]を計算する。

関数 solverNewmarkBeta とそれを用いた数値計算の簡単な説明

    an1_num = fn1 - c*(vn + dt*an/2.0) - k*(xn + dt*vn + (1.0/2.0 - beta)*(dt**2)*an)
    an1_den = m + dt*c/2.0 + beta*(dt**2)*k
    an1 = an1_num/an1_den
    vn1 = vn + dt*(an + an1)/2.0
    xn1 = xn + dt*vn + (dt**2)*an/2.0 + beta*(dt**2)*(an1 - an)

こちらでは、ニューマークβ法の直接計算法にて示した具体的計算(1自由度の場合)の結果で得られた計算式をそのまま実装している。変数の意味は関数definiteSolverNewmarkBetaにて示したとおりである。
こちらの式を使用した数値計算部分は以下の部分である。

# ニューマークβ法による計算
a_nm[i],v_nm[i], x_nm[i] = solverNewmarkBeta(m,c,k,
                                             x_nm[i-1],
                                             v_nm[i-1],
                                             a_nm[i-1],
                                             f,dt,b)

プログラムの実行結果

実行結果を以下に示す。

解析解を黒実線、シンボリック計算にて算出した式を利用したdefiniteSolverNewmarkBetaを使用した数値計算結果は緑の実線、直接計算式を実装したsolverNewmarkBeta数値計算結果を赤点線で示している。実装が間違っていなければ両者は一致するはずであるが、一致しているため正しく実装できていると思われる。また、解析解とも概ね一致している。

参考文献

戸川隼人著.「有限要素法による振動解析」サイエンス社.1975年