Newmark-β 法で1次元の波動方程式

微分の表記

微分の表記を短く書けるものに変えます。時間の微分にニュートンの表記を使います。

u˙=ut,u¨=2ut2, \dot{u} = \frac{\partial u}{\partial t} ,\quad \ddot{u} = \frac{\partial^2 u}{\partial t^2} ,\quad \ldots

空間方向の微分にラグランジュの表記を使います。

u=ux,u=2ux2, u' = \frac{\partial u}{\partial x} ,\quad u'' = \frac{\partial^2 u}{\partial x^2} ,\quad \ldots

離散化で出てくる差分を Δt\Delta_tΔx\Delta_x と表記します。以前の記事では dtdtdxdx と書いていましたが、乗算が abcdta b c dt のようになると区切りが分からなくなるので変えました。

uu のインデックスを下付き文字で表記します。

un=u(x,t+nΔt)un+1=u(x,t+(n+1)Δt)un,x+1=u(x+Δx,t+nΔt)un,x1=u(xΔx,t+(n+1)Δt) \begin{aligned} u_{n} &= u(x, t + n\,\Delta_t)\\ u_{n+1} &= u(x, t + (n + 1)\,\Delta_t)\\ u_{n, x+1} &= u(x + \Delta_x, t + n\,\Delta_t)\\ u_{n, x-1} &= u(x - \Delta_x, t + (n + 1)\,\Delta_t)\\ \end{aligned}

波動方程式

1次元の波動方程式です。

u¨=c2u \ddot{u} = c^2 u''

波動方程式はエネルギーの減衰が考慮されていないので時間が経っても波が減衰しないはずですが Implicit FDM で離散化すると減衰が起こります。 Schweickart, James, Marschner の "Animating Elastic Rods with Sound" に習って Newmark-β 法を使うことでエネルギーの減衰を改善します。

Newmark-β 法

Newmark-β 法では速度と位置の計算に次の式を使います。

u˙n+1=u˙n+Δt((1γ)u¨n+γu¨n+1)un+1=un+Δtu˙n+Δt2((12β)u¨n+βu¨n+1) \begin{aligned} \dot{u}_{n+1} &= \dot{u}_{n} + \Delta_t \left( (1 - \gamma) \ddot{u}_{n} + \gamma \ddot{u}_{n+1} \right)\\ u_{n+1} &= u_{n} + \Delta_t \dot{u}_{n} + \Delta_t^2 \left( \left( \frac{1}{2} - \beta \right) \ddot{u}_{n} + \beta \ddot{u}_{n+1} \right) \end{aligned}

β\betaγ\gamma が調整できる形になっていますが γ\gamma については 1/21/2 以外の値にすると精度が下がるそうです。 β\beta の値はWikipediaの記事1/41/4 、 Gavin の "Numerical Integration in Structural Dynamics"1/61/6 の場合が紹介されていました。この文章では γ\gamma1/21/2 で固定して β\beta は調整できる形にしています。

Newmark-β 法による波動方程式の離散化

位置 uu 、速度 u˙\dot{u} 、加速度 u¨\ddot{u} からなる連立微分方程式を立てます。位置と速度の式は Newmark-β 法のもので、加速度の式は1次元の波動方程式です。

u¨n+1=c2un+1u˙n+1=u˙n+Δt2(u¨n+u¨n+1)un+1=un+Δtu˙n+Δt2((12β)u¨n+βu¨n+1) \begin{aligned} \ddot{u}_{n+1} &= c^2 u''_{n+1}\\ \dot{u}_{n+1} &= \dot{u}_{n} + \frac{\Delta_t}{2} \left( \ddot{u}_{n} + \ddot{u}_{n+1} \right)\\ u_{n+1} &= u_{n} + \Delta_t \dot{u}_{n} + \Delta_t^2 \left( \left( \frac{1}{2} - \beta \right) \ddot{u}_{n} + \beta \ddot{u}_{n+1} \right) \end{aligned}

u¨n+1\ddot{u}_{n+1} の右辺について un+1u''_{n+1} を有限差分に変形して un+1u_{n+1} の右辺を代入します。

u¨n+1=c2un+1=c2Δx2(un+1,x12un+1+un+1,x+1)=c2Δx2((un,x12un+un,x+1)+Δt(u˙n,x12u˙n+u˙n,x+1)+Δt2(12β)(u¨n,x12u¨n+u¨n,x+1)+Δt2β(u¨n+1,x12u¨n+1+u¨n+1,x+1)) \begin{aligned} \ddot{u}_{n+1} &= c^2 u''_{n+1}\\ &= \frac{c^2}{\Delta_x^2} (u_{n+1, x-1} - 2 u_{n+1} + u_{n+1, x+1})\\ &= \frac{c^2}{\Delta_x^2} \Biggl( (u_{n, x-1} - 2 u_{n} + u_{n, x+1}) + \Delta_t ( \dot{u}_{n, x-1} - 2 \dot{u}_{n} + \dot{u}_{n, x+1})\\ &\qquad\qquad+ \Delta_t^2 \left( \frac{1}{2} - \beta \right) ( \ddot{u}_{n, x-1} - 2 \ddot{u}_{n} + \ddot{u}_{n, x+1})\\ &\qquad\qquad+ \Delta_t^2 \beta ( \ddot{u}_{n+1, x-1} - 2 \ddot{u}_{n+1} + \ddot{u}_{n+1, x+1}) \Biggr)\\ \end{aligned}

n+1n+1 の項を左辺に移項します。

(Δx2c2+2Δt2β)u¨n+1Δt2β(u¨n+1,x1+u¨n+1,x+1)=(un,x12un+un,x+1)+Δt(u˙n,x12u˙n+u˙n,x+1)+Δt2(12β)(u¨n,x12u¨n+u¨n,x+1) \begin{aligned} & \left(\frac{\Delta_x^2}{c^2} + 2 \Delta_t^2 \beta \right) \ddot{u}_{n+1} - \Delta_t^2 \beta ( \ddot{u}_{n+1, x-1} + \ddot{u}_{n+1, x+1})\\ &= (u_{n, x-1} - 2 u_{n} + u_{n, x+1}) + \Delta_t ( \dot{u}_{n, x-1} - 2 \dot{u}_{n} + \dot{u}_{n, x+1})\\ &\quad+ \Delta_t^2 \left( \frac{1}{2} - \beta \right) ( \ddot{u}_{n, x-1} - 2 \ddot{u}_{n} + \ddot{u}_{n, x+1})\\ \end{aligned}

整理します。

C0u¨n+1+C1(u¨n+1,x1+u¨n+1,x+1)=(un,x12un+un,x+1)+Δt(u˙n,x12u˙n+u˙n,x+1)+C2(u¨n,x12u¨n+u¨n,x+1) \begin{aligned} C_0 \ddot{u}_{n+1} + C_1 ( \ddot{u}_{n+1, x-1} + \ddot{u}_{n+1, x+1}) &= (u_{n, x-1} - 2 u_{n} + u_{n, x+1})\\ &\quad+ \Delta_t ( \dot{u}_{n, x-1} - 2 \dot{u}_{n} + \dot{u}_{n, x+1})\\ &\quad+ C_2 ( \ddot{u}_{n, x-1} - 2 \ddot{u}_{n} + \ddot{u}_{n, x+1})\\ \end{aligned}

境界となる左右の端では式の形が変わります。次の式は左端の式です。 x+1x+1x1x-1 に変えると右端の式になります。 b=1b=1 のとき固定端、 b=2b=2 のとき自由端になります。

C0u¨n+1+bC1u¨n+1,x+1=(bun,x+12un)+Δt(bu˙n,x+12u˙n)+C2(bu¨n,x+12u¨n) \begin{aligned} C_0 \ddot{u}_{n+1} + b\, C_1 \ddot{u}_{n+1, x+1} &= (b\, u_{n, x+1} - 2 u_{n})\\ &\quad+ \Delta_t (b\, \dot{u}_{n, x+1} - 2 \dot{u}_{n})\\ &\quad+ C_2 (b\, \ddot{u}_{n, x+1} - 2 \ddot{u}_{n}) \end{aligned}

連立方程式を立てます。

{Au¨n+1=un+Δtu˙n+C2u¨nun+1=un+Δtu˙n+C2u¨nC1u¨n+1u˙n+1=u˙n+C3(u¨n+u¨n+1),A=[C0bLC10000C1C0C10000C1C0C1000000bRC1C0]C0=Δx2c2+2Δt2β,C1=Δt2β,C2=Δt22+C1,C3=Δt2 \begin{aligned} \begin{cases} \mathbf{A} \ddot{\mathbf{u}}_{n+1} = \mathbf{u}''_n + \Delta_t \dot{\mathbf{u}}''_{n} + C_2 \ddot{\mathbf{u}}''_{n}\\ u_{n+1} = u_{n} + \Delta_t \dot{u}_{n} + C_2 \ddot{u}_{n} - C_1 \ddot{u}_{n+1}\\ \dot{u}_{n+1} = \dot{u}_{n} + C_3 \left( \ddot{u}_{n} + \ddot{u}_{n+1} \right)\\ \end{cases} \;,&\quad \mathbf{A} = \begin{bmatrix} C_0 & b_L C_1 & 0 & 0 & \cdots & 0 & 0\\ C_1 & C_0 & C_1 & 0 & \cdots & 0 & 0\\ 0 & C_1 & C_0 & C_1 & \cdots & 0 & 0\\ & \vdots & & & & \vdots &\\ 0 & 0 & 0 & 0 & \cdots & b_R C_1 & C_0 \end{bmatrix}\\ C_0 = \frac{\Delta_x^2}{c^2} + 2 \Delta_t^2 \beta ,\quad C_1 = - \Delta_t^2 \beta ,&\quad C_2 = \frac{\Delta_t^2}{2} + C_1 ,\quad C_3 = \frac{\Delta_t}{2} \end{aligned}

実装はGitHubの別ページに分けました。リンク先の Wave1DNewmarkBeta クラスになります。

実装を読む (github.com)

Newmark-β 法によるばね-ダンパ波動方程式の離散化

波動方程式にばねの項 kuk u とダンパの項 au˙a \dot{u} を加えることで波の減衰を調整できるようにします。この文章では次の式をばね-ダンパ波動方程式と呼ぶことにします。

u¨+au˙+ku=c2u \ddot{u} + a \dot{u} + k u = c^2 u''

Newmark-β 法の連立方程式を立てます。

u¨n+1=c2un+1au˙n+1+kun+1u˙n+1=u˙n+Δt2(u¨n+u¨n+1)un+1=un+Δtu˙n+Δt2((12β)u¨n+βu¨n+1) \begin{aligned} \ddot{u}_{n+1} &= c^2 u''_{n+1} - a \dot{u}_{n+1} + k u_{n+1}\\ \dot{u}_{n+1} &= \dot{u}_{n} + \frac{\Delta_t}{2} \left( \ddot{u}_{n} + \ddot{u}_{n+1} \right)\\ u_{n+1} &= u_{n} + \Delta_t \dot{u}_{n} + \Delta_t^2 \left( \left( \frac{1}{2} - \beta \right) \ddot{u}_{n} + \beta \ddot{u}_{n+1} \right) \end{aligned}

u¨n+1\ddot{u}_{n+1} の右辺について u˙n+1\dot{u}_{n+1}un+1u_{n+1} を代入します。

u¨n+1=c2Δx2((un,x12un+un,x+1)+Δt(u˙n,x12u˙n+u˙n,x+1)+Δt2(12β)(u¨n,x12u¨n+u¨n,x+1)+Δt2β(u¨n+1,x12u¨n+1+u¨n+1,x+1))a(u˙n+Δt2(u¨n+u¨n+1))k(un+Δtu˙n+Δt2(12β)u¨n+Δt2βu¨n+1) \begin{aligned} \ddot{u}_{n+1} &= \frac{c^2}{\Delta_x^2} \Biggl( (u_{n, x-1} - 2 u_{n} + u_{n, x+1}) + \Delta_t ( \dot{u}_{n, x-1} - 2 \dot{u}_{n} + \dot{u}_{n, x+1})\\ &\qquad\qquad+ \Delta_t^2 \left( \frac{1}{2} - \beta \right) ( \ddot{u}_{n, x-1} - 2 \ddot{u}_{n} + \ddot{u}_{n, x+1})\\ &\qquad\qquad+ \Delta_t^2 \beta ( \ddot{u}_{n+1, x-1} - 2 \ddot{u}_{n+1} + \ddot{u}_{n+1, x+1}) \Biggr)\\ &\quad- a \Biggl( \dot{u}_{n} + \frac{\Delta_t}{2} \left( \ddot{u}_{n} + \ddot{u}_{n+1} \right) \Biggr)\\ &\quad- k \Biggl( u_{n} + \Delta_t \dot{u}_{n} + \Delta_t^2 \left( \frac{1}{2} - \beta \right) \ddot{u}_{n} + \Delta_t^2 \beta \ddot{u}_{n+1} \Biggr) \end{aligned}

n+1n+1 の項を左辺に移項します。

(1+aΔt2+kΔt2β+2c2Δx2Δt2β)u¨n+1c2Δx2Δt2β(u¨n+1,x1+u¨n+1,x+1)=c2Δx2((un,x1+un,x+1)+Δt(u˙n,x1+u˙n,x+1)+Δt2(12β)(u¨n,x1+u¨n,x+1))(aΔt2+kΔt2(12β)+2c2Δx2Δt2(12β))u¨n(a+kΔt+2c2Δx2Δt)u˙n(k+2c2Δx2)un \begin{aligned} & ( 1 + a \frac{\Delta_t}{2} + k \Delta_t^2 \beta + 2 \frac{c^2}{\Delta_x^2} \Delta_t^2 \beta ) \ddot{u}_{n+1} - \frac{c^2}{\Delta_x^2} \Delta_t^2 \beta ( \ddot{u}_{n+1, x-1} + \ddot{u}_{n+1, x+1})\\ &= \frac{c^2}{\Delta_x^2} \Biggl( (u_{n, x-1} + u_{n, x+1}) + \Delta_t (\dot{u}_{n, x-1} + \dot{u}_{n, x+1}) + \Delta_t^2 \left( \frac{1}{2} - \beta \right) ( \ddot{u}_{n, x-1} + \ddot{u}_{n, x+1}) \Biggr)\\ &\quad- \Biggl( a \frac{\Delta_t}{2} + k \Delta_t^2 \left( \frac{1}{2} - \beta \right) + 2 \frac{c^2}{\Delta_x^2} \Delta_t^2 \left( \frac{1}{2} - \beta \right) \Biggr) \ddot{u}_{n}\\ &\quad- \Biggl( a + k \Delta_t + 2 \frac{c^2}{\Delta_x^2} \Delta_t \Biggr) \dot{u}_{n} \\ &\quad- \Biggl( k + 2 \frac{c^2}{\Delta_x^2} \Biggr) u_{n} \end{aligned}

整理します。

C0u¨n+1+C1(u¨n+1,x1+u¨n+1,x+1)=C2((un,x1+un,x+1)+Δt(u˙n,x1+u˙n,x+1)+C3(u¨n,x1+u¨n,x+1))C4u¨nC5u˙nC6un \begin{aligned} & C_0 \ddot{u}_{n+1} + C_1 ( \ddot{u}_{n+1, x-1} + \ddot{u}_{n+1, x+1})\\ &= C_2 \bigl( (u_{n, x-1} + u_{n, x+1}) + \Delta_t (\dot{u}_{n, x-1} + \dot{u}_{n, x+1}) + C_3 (\ddot{u}_{n, x-1} + \ddot{u}_{n, x+1}) \bigr)\\ &\quad- C_4 \ddot{u}_{n} - C_5 \dot{u}_{n} - C_6 u_{n}\\ \end{aligned}

境界となる左右の端では式の形が変わります。次の式は左端の式です。 x+1x+1x1x-1 に変えると右端の式になります。 b=1b=1 のとき固定端、 b=2b=2 のとき自由端になります。

C0u¨n+1+bC1u¨n+1,x+1=bC2(un,x+1+Δtu˙n,x+1+C3u¨n,x+1)C4u¨nC5u˙nC6un \begin{aligned} C_0 \ddot{u}_{n+1} + b\,C_1 \ddot{u}_{n+1, x+1} &= b\,C_2 \bigl( u_{n, x+1} + \Delta_t \dot{u}_{n, x+1} + C_3 \ddot{u}_{n, x+1} \bigr)\\ &\quad- C_4 \ddot{u}_{n} - C_5 \dot{u}_{n} - C_6 u_{n}\\ \end{aligned}

連立方程式を立てます。

C0=1+aC7+C6C8C1=C2C8C2=c2/Δx2C3=Δt2(1/2β)C4=C3C6+aC7C5=a+ΔtC6C6=k+2C2C7=Δt/2C8=Δt2β{Au¨n+1=bu˙n+1=u˙n+C7(u¨n+u¨n+1)un+1=un+Δtu˙n+C3u¨n+C8u¨n+1A=[C0bLC10000C1C0C10000C1C0C1000000bRC1C0]b=C2((un,x1+un,x+1)+Δt(u˙n,x1+u˙n,x+1)+C3(u¨n,x1+u¨n,x+1))C4u¨nC5u˙nC6un \begin{aligned} C_0 &= 1 + a C_7 + C_6 C_8\\ C_1 &= - C_2 C_8\\ C_2 &= c^2 / \Delta_x^2\\ C_3 &= \Delta_t^2 \left( 1 / 2 - \beta \right)\\ C_4 &= C_3 C_6 + a C_7\\ C_5 &= a + \Delta_t C_6\\ C_6 &= k + 2 C_2\\ C_7 &= \Delta_t / 2\\ C_8 &= \Delta_t^2 \beta\\ \end{aligned} \qquad \begin{aligned} &\quad \begin{cases} \mathbf{A} \ddot{\mathbf{u}}_{n+1} = \mathbf{b}\\ \dot{u}_{n+1} = \dot{u}_{n} + C_7 \left( \ddot{u}_{n} + \ddot{u}_{n+1} \right)\\ u_{n+1} = u_{n} + \Delta_t \dot{u}_{n} + C_3 \ddot{u}_{n} + C_8 \ddot{u}_{n+1}\\ \end{cases}\\ \\ \mathbf{A} &= \begin{bmatrix} C_0 & b_L C_1 & 0 & 0 & \cdots & 0 & 0\\ C_1 & C_0 & C_1 & 0 & \cdots & 0 & 0\\ 0 & C_1 & C_0 & C_1 & \cdots & 0 & 0\\ & \vdots & & & & \vdots &\\ 0 & 0 & 0 & 0 & \cdots & b_R C_1 & C_0 \end{bmatrix}\\ \\ \mathbf{b} &= C_2 \bigl( (\mathbf{u}_{n, x-1} + \mathbf{u}_{n, x+1}) + \Delta_t (\dot{\mathbf{u}}_{n, x-1} + \dot{\mathbf{u}}_{n, x+1})\\ &\qquad+ C_3 (\ddot{\mathbf{u}}_{n, x-1} + \ddot{\mathbf{u}}_{n, x+1}) \bigr) - C_4 \ddot{\mathbf{u}}_{n} - C_5 \dot{\mathbf{u}}_{n} - C_6 \mathbf{u}_{n}\\ \end{aligned}

実装はGitHubの別ページに分けました。リンク先の DampedWave1DNewmarkBeta クラスになります。

実装を読む (github.com)

デモ

キャンバスをクリックすると波が起きます。

Δt=1/60\Delta_t = 1/60 で固定しています。

ソルバにガウス-ザイデル法を使っています。 cΔt/Δxc \Delta_t / \Delta_x が大きいとソルバの収束が遅くなって発散します。ガウス-ザイデル法の反復を増やすことで発散を抑えています。

β<1/4\beta < 1/4 のとき、他のパラメータの値によっては発散します。

参考サイト