Newmark Beta法
参考：
微分方程式の数値計算 線形加速度法(ニューマークβ法)
台形公式(Wikipedia)
Newmark-beta method(Wikipedia)
平均値の定理・ロルの定理とその証明
構造物を動的解析するために、構造物に力$F$がかかった時、$\Delta t$秒後の構造物の変位$x$を知りたい
運動方程式$m\ddot{x}=F$から、加速度は$\ddot{x}=\frac{F}{m}$で求められる
さらに、加速度$\ddot{x}$を積分すると、速度$\dot{x}$
速度$\dot{x}$を積分すると変位$x$
これで変位$x$が求められるが、プログラムは積分ができないので、
数値積分という近似によってその解を求める
Newmarkβ法は、その一種である
§I - 加速度の数値積分
等加速度直線運動の公式を考える（加速度の不定積分($C=\dot{x}_0$)）
$$
\begin{aligned}
\dot{x}&=\int \ddot{x}dt=\ddot{x}t+C\\
\dot{x}&=\dot{x}_0+\ddot{x}t
\end{aligned}
$$
これは加速度が水平の話で、$\ddot{x}$は定数
加速度とは、つまり単位時間における速度の増加量であるため、時間$t$を掛けることで、その時間での速度の増加量がわかる。
これが$\ddot{x}t$の表すところであり、加速度曲線(下図)の面積にあたる。
$\dot{x}$は、 計算前の速度に時間での速度の増加量が足すことで求められる。
台形則
加速度が時間変化する場合、$\ddot{x}$は定数ではなく、$t$の値で変化する方程式$\ddot{x}(t)$になる
数値積分では、非線形方程式に対応するため、曲線を直線の集まりで近似してコンピュータが計算できる形にする(下図)
厳密には、$t_b\sim (t_b+\Delta t)$の間では線形変化していくものと仮定する
当然、$\Delta t$が小さくなるほど、誤差は減少する
$\Delta t$秒後の速度$\dot{x}(t_b+\Delta t)$は、現在の速度$\dot{x}(t_b)$と加速度の増加量$\dot{x}(\Delta t)$の和である
現在の速度$\dot{x}(t_b)$は過去の計算の総和なので、既知として、
速度の増加量$\dot{x}(\Delta t)$は、先の図の$\ddot{x}t$の部分にあたり、加速度曲線における面積は台形の公式（$\frac{上底+下底}{2}\cdot高さ$）を使って求められる
これを台形則と呼んだりする
$$
\begin{aligned}
\dot{x}(t_b+\Delta t)&=\dot{x}(t_b)+\dot{x}(\Delta t)\\
\dot{x}(t_b+\Delta t)&=\dot{x}(t_b)+\frac{\ddot{x}(t_b)+\ddot{x}(t_b+\Delta t)}{2}\Delta t
\end{aligned}
$$
物理的な観点においてはこれでよいのだが、流体解析では、計算安定性を制御するために、数値粘性項を定義することがある
Newmarkβ法の利点は数値粘性項を定義することで式の性質を簡単に変更できる点にある（たぶん）
流体解析的な視点で、計算安定性を確保するために人為的に与えられる値のこと。流体の移流における粘性に新たな項を加えていることに相当する
内分点の定義
加速度曲線において、$\ddot{x}(t_b)$と$\ddot{x}(t_b+\Delta t)$を内分する$\ddot{x}(\gamma)$を考える
内分する点は下図のような公式で求める
$\frac{n}{m+n},\frac{m}{m+n}$はつまるところA-B間をPで分割したときの割合である
よって下図の公式は1変数で以下のように書ける
$$
\begin{aligned}
P=(1-\gamma)A+\gamma B\\
where\space 0<\gamma<1
\end{aligned}
$$
引用元︰Try IT
加速度曲線の点に書き直すと、
$$
\begin{aligned}
&\ddot{x}(\gamma)=(1-\gamma)\ddot{x}(t)+\gamma\ddot{x}(t+\Delta t)\\
&where\space 0<\gamma<1
\end{aligned}
$$
$\gamma$を小さくするほど$\ddot{x}(t)$に近づき、粘性が強くなる
(加速度が爆発したときに安定性を保持できるが、計算精度が下がる)
$\gamma=0.5$とすれば、台形公式の$\frac{\ddot{x}(t)+\ddot{x}(t+\Delta t)}{2}$と同じで、ほとんどの計算手法で$\gamma$は$0.5$である
ということで、速度は次式のように表される
$$
\begin{aligned}
\dot{x}(t_b+\Delta t)&=\dot{x}(t_b)+\ddot{x}(\gamma)\Delta t\\
\dot{x}(t_b+\Delta t)&=\dot{x}(t_b)+((1-\gamma)\ddot{x}(t_b)+\gamma\ddot{x}(t_b+\Delta t))\Delta t\\
&where\space 0<\gamma<1
\end{aligned}
$$§II - 変位の数値積分
§Iと同じ方法で解くと、平均加速度法と呼ばれる別の手法になるが、とりあえず解いてみる
速度のときと同じように、台形則で変位を求めると、次式のようになる
（加速度→速度、速度→変位になっただけ）
$$
\begin{aligned}
x(t_b+\Delta t)&=x(t_b)+x(\Delta t)\\
x(t_b+\Delta t)&=x(t_b)+\frac{\dot{x}(t_b)+\dot{x}(t_b+\Delta t)}{2}\Delta t
\end{aligned}
$$
さらに、前のセクションの台形則で求めた$\dot{x}(t_b+\Delta t)$を代入すると、
$$
x(t_b+\Delta t)=x(t_b)+\frac{\dot{x}(t_b)+\left(\dot{x}(t_b)+\frac{\ddot{x}(t_b)+\ddot{x}(t_b+\Delta t)}{2}\Delta t\right)}{2}\Delta t
$$
整理して、
$$
x(t_b+\Delta t)=x(t_b)+\dot{x}(t_b)\Delta t+\frac{1}{2}\cdot\frac{\ddot{x}(t_b)+\ddot{x}(t_b+\Delta t)}{2}(\Delta t)^2
$$βを定義するには、線形加速度法の解き方を用いる
具体的には、テイラー展開を使う
テイラー展開は関数$f$において、$a$から少し離れた$x$を求める方法で、項が多ければ多いほどより遠くまで求められる
$$
\sum_{n=0}^{\infty}\frac{f^{(n)}(a)}{n!}(x-a)^n
$$
$x(t_b+\Delta t)$について、3次までテイラー展開すると、
$a=t_b$
$x=t_b+\Delta t$
$$
x(t_b+\Delta t)=x(t_b)+\frac{\dot{x}(t_b)}{1!}\Delta t+\frac{\ddot{x}(t_b)}{2!}(\Delta t)^2+\frac{\dddot{x}(t_b)}{3!}(\Delta t)^3
$$
ここで、数値積分（というより線形加速度法）では、加速度の傾き$\dddot{x}(t_b)$は$\Delta t$の中で、線形変化として扱うので、
$$
\dddot{x}(t_b)=\frac{\ddot{x}(t_b+\Delta t)-\ddot{x}(t_b)}{\Delta t}
$$
これを代入して、
$$
x(t_b+\Delta t)=x(t_b)+\frac{\dot{x}(t_b)}{1!}\Delta t+\frac{\ddot{x}(t_b)}{2!}(\Delta t)^2+\frac{1}{3!}\cdot\frac{\ddot{x}(t_b+\Delta t)-\ddot{x}(t_b)}{\Delta t}(\Delta t)^3
$$
線形加速度法では、3次項の係数が$\frac{1}{3!}$である
この3次項の係数を変えると式は様々な性質に変化する
そのため、$\frac{1}{3!}$を$\beta$として変更できるようにしたのがNewmarkβ法
$$
x(t_b+\Delta t)=x(t_b)+\frac{\dot{x}(t_b)}{1!}\Delta t+\frac{\ddot{x}(t_b)}{2!}(\Delta t)^2+\beta\frac{\ddot{x}(t_b+\Delta t)-\ddot{x}(t_b)}{\Delta t}(\Delta t)^3
$$
これを整理
分母を$\Delta t$で消去し、$\beta,(\Delta t)^2$を分配
$$
x(t_b+\Delta t)=x(t_b)+\frac{\dot{x}(t_b)}{1!}\Delta t+\frac{\ddot{x}(t_b)}{2!}(\Delta t)^2+\beta\ddot{x}(t_b+\Delta t)(\Delta t)^2-\beta\ddot{x}(t_b)(\Delta t)^2
$$
全体に2をかける
$$
2x(t_b+\Delta t)=2x(t_b)+2\dot{x}(t_b)\Delta t+\ddot{x}(t_b)(\Delta t)^2+2\beta\ddot{x}(t_b+\Delta t)(\Delta t)^2-2\beta\ddot{x}(t_b)(\Delta t)^2
$$
$\ddot{x}(t_b)(\Delta t)^2$でまとめる
$$
2x(t_b+\Delta t)=2x(t_b)+2\dot{x}(t_b)\Delta t+(1-2\beta)\ddot{x}(t_b)(\Delta t)^2+2\beta\ddot{x}(t_b+\Delta t)(\Delta t)^2
$$
$(\Delta t)^2$でまとめる
$$
2x(t_b+\Delta t)=2x(t_b)+2\dot{x}(t_b)\Delta t+((1-2\beta)\ddot{x}(t_b)+2\beta\ddot{x}(t_b+\Delta t))(\Delta t)^2
$$
全体を2で割る
$$
\begin{aligned}
x(t_b+\Delta t)&=x(t_b)+\dot{x}(t_b)\Delta t+\frac{1}{2}((1-2\beta)\ddot{x}(t_b)+2\beta\ddot{x}(t_b+\Delta t))(\Delta t)^2\\
&where\space 0<\beta<\frac{1}{2}
\end{aligned}
$$
βを変えることで式の性質が変わる
$\beta=0$：陽解法
$\beta=\frac{1}{6}$：線形加速度法
$\beta=\frac{1}{4}$：平均加速度法
$\beta=\frac{1}{2}$：陰解法
§III - まとめ
Newmarkβ法では、変位、速度を次式のように表す
$$
\begin{aligned}
x(t_b+\Delta t)&=x(t_b)+\dot{x}(t_b)\Delta t+\frac{1}{2}((1-2\beta)\ddot{x}(t_b)+2\beta\ddot{x}(t_b+\Delta t))(\Delta t)^2\\
&where\space 0<\beta<\frac{1}{2}\\
\\
\dot{x}(t_b+\Delta t)&=\dot{x}(t_b)+((1-\gamma)\ddot{x}(t_b)+\gamma\ddot{x}(t_b+\Delta t))\Delta t\\
&where\space 0<\gamma<1
\end{aligned}
$$速度は、台形則で求めたが、もちろん変位と同じようにテイラー展開で求められる
$\dot{x}(t_b+\Delta t)$について、3次までテイラー展開すると、
$$
\dot{x}(t_b+\Delta t)=\dot{x}(t_b)+\frac{\ddot{x}(t_b)}{1!}\Delta t+\frac{\dddot{x}(t_b)}{2!}(\Delta t)^2
$$
例のごとく$\dddot{x}(t_b)=\frac{\ddot{x}(t_b+\Delta t)-\ddot{x}(t_b)}{\Delta t}$とすると、
$$
\begin{aligned}
\dot{x}(t_b+\Delta t)&=\dot{x}(t_b)+\ddot{x}(t_b)\Delta t+\frac{1}{2}\left(\ddot{x}(t_b+\Delta t)-\ddot{x}(t_b)\right)\Delta t\\
\dot{x}(t_b+\Delta t)&=\dot{x}(t_b)+\ddot{x}(t_b)\Delta t+\frac{1}{2}\ddot{x}(t_b+\Delta t)\Delta t-\frac{1}{2}\ddot{x}(t_b)\Delta t\\
\dot{x}(t_b+\Delta t)&=\dot{x}(t_b)+(1-\frac{1}{2})\ddot{x}(t_b)\Delta t+\frac{1}{2}\ddot{x}(t_b+\Delta t)\Delta t\\
\dot{x}(t_b+\Delta t)&=\dot{x}(t_b)+\left((1-\frac{1}{2})\ddot{x}(t_b)+\frac{1}{2}\ddot{x}(t_b+\Delta t)\right)\Delta t
\end{aligned}
$$
要は、この$\frac{1}{2}$が$\gamma$である