2次元のばね-ダンパ波動方程式

2次元のばね-ダンパ波動方程式を実装します。

表記について

時間のインデックスの表記を nn から tt に変えます。 nn は次元の数に使うことにします。

ラプラシアンの離散化

nn 次元のラプラシアン有限差分の形で離散化します。

2ut=1Δx2(neighbor(ut)2nut)neighbor(ut)=iaxis(n)(ut,i1+ut,i+1),axis(n)={x,y,z,} \begin{aligned} \nabla^{2} u_{t} &= \frac{1}{\Delta_x^{2}} \left( \mathtt{neighbor}(u_{t}) - 2n\,u_{t} \right)\\ \mathtt{neighbor}(u_{t}) &= \sum_{i\,\in \mathtt{axis}(n)} (u_{t,i-1} + u_{t,i+1}) ,\qquad \mathtt{axis}(n) = \{x, y, z, \dots \} \end{aligned}

neighbor(ut)\mathtt{neighbor}(u_{t})近傍の総和です。 axis(n)\mathtt{axis}(n)nn 次元の各軸の集合を表しているつもりですが、フォーマルな書き方がよく分からないのでごまかして書いています。

離散化

後で次元を変えられる形で離散化します。 nn 次元のばね-ダンパ波動方程式です。

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

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

u¨t+1=c22ut+1au˙t+1+kut+1u˙t+1=u˙t+Δt2(u¨t+u¨t+1)ut+1=ut+Δtu˙t+Δt2((12β)u¨t+βu¨t+1) \begin{aligned} \ddot{u}_{t+1} &= c^{2} \nabla^{2} u_{t+1} - a \dot{u}_{t+1} + k u_{t+1}\\ \dot{u}_{t+1} &= \dot{u}_{t} + \frac{\Delta_t}{2} \left( \ddot{u}_{t} + \ddot{u}_{t+1} \right)\\ u_{t+1} &= u_{t} + \Delta_t \dot{u}_{t} + \Delta_t^{2} \left( \left( \frac{1}{2} - \beta \right) \ddot{u}_{t} + \beta \ddot{u}_{t+1} \right) \end{aligned}

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

u¨t+1=c22(ut+Δtu˙t+Δt2((12β)u¨t+βu¨t+1))au˙t+1+kut+1=c2Δx2((neighbor(ut)2nut)+Δt(neighbor(u˙t)2nu˙t)+Δt2(12β)(neighbor(u¨t)2nu¨t)+Δt2β(neighbor(u¨t+1)2nu¨t+1))a(u˙t+Δt2(u¨t+u¨t+1))k(ut+Δtu˙t+Δt2(12β)u¨t+Δt2βu¨t+1) \begin{aligned} \ddot{u}_{t+1} &= c^{2} \nabla^{2} \left( u_{t} + \Delta_t \dot{u}_{t} + \Delta_t^{2} \left( \left( \frac{1}{2} - \beta \right) \ddot{u}_{t} + \beta \ddot{u}_{t+1} \right) \right) - a \dot{u}_{t+1} + k u_{t+1}\\ &= \frac{c^{2}}{\Delta_x^{2}} \Biggl( ( \mathtt{neighbor}(u_{t}) - 2n\,u_{t}) + \Delta_t ( \mathtt{neighbor}(\dot{u}_{t}) - 2n\,\dot{u}_{t})\\ &\qquad\qquad+ \Delta_t^{2} \left( \frac{1}{2} - \beta \right) ( \mathtt{neighbor}(\ddot{u}_{t}) - 2n\,\ddot{u}_{t})\\ &\qquad\qquad+ \Delta_t^{2} \beta ( \mathtt{neighbor}(\ddot{u}_{t+1}) - 2n\,\ddot{u}_{t+1}) \Biggr)\\ &\quad- a \Biggl( \dot{u}_{t} + \frac{\Delta_t}{2} \left( \ddot{u}_{t} + \ddot{u}_{t+1} \right) \Biggr)\\ &\quad- k \Biggl( u_{t} + \Delta_t \dot{u}_{t} + \Delta_t^{2} \left( \frac{1}{2} - \beta \right) \ddot{u}_{t} + \Delta_t^{2} \beta \ddot{u}_{t+1} \Biggr) \end{aligned}

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

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

整理します。

C0u¨t+1+C1neighbor(u¨t+1)=C2(neighbor(ut)+Δtneighbor(u˙t)+C3neighbor(u¨t))C4u¨tC5u˙tC6ut \begin{aligned} & C_0 \ddot{u}_{t+1} + C_1 \,\mathtt{neighbor}(\ddot{u}_{t+1})\\ &= C_2 \bigl( \mathtt{neighbor}(u_{t}) + \Delta_t \,\mathtt{neighbor}(\dot{u}_{t}) + C_3 \,\mathtt{neighbor}(\ddot{u}_{t}) \bigr)\\ &\quad- C_4 \ddot{u}_{t} - C_5 \dot{u}_{t} - C_6 u_{t}\\ \end{aligned}

連立方程式を立てます。1次元のときと比べると、定数については C6C_6 だけが変わっています。

C0=1+aC7+C6C8C1=C2C8C2=c2/Δx2C3=Δt2(1/2β)C4=C3C6+aC7C5=a+ΔtC6C6=k+2nC2C7=Δt/2C8=Δt2β{Au¨t+1=bu˙t+1=u˙t+C7(u¨t+u¨t+1)ut+1=ut+Δtu˙t+C3u¨t+C8u¨t+1b=C2(neighbor(ut)+Δtneighbor(u˙t)+C3neighbor(u¨t))C4u¨tC5u˙tC6ut \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 + 2n\,C_2\\ C_7 &= \Delta_t / 2\\ C_8 &= \Delta_t^{2} \beta\\ \end{aligned} \quad \begin{aligned} &\quad \begin{cases} \mathbf{A} \ddot{\mathbf{u}}_{t+1} = \mathbf{b}\\ \dot{u}_{t+1} = \dot{u}_{t} + C_7 \left( \ddot{u}_{t} + \ddot{u}_{t+1} \right)\\ u_{t+1} = u_{t} + \Delta_t \dot{u}_{t} + C_3 \ddot{u}_{t} + C_8 \ddot{u}_{t+1}\\ \end{cases}\\ \\ \mathbf{b} &= C_2 \bigl( \mathtt{neighbor}(\mathbf{u}_{t}) + \Delta_t \mathtt{neighbor}(\dot{\mathbf{u}}_{t}) + C_3 \mathtt{neighbor}(\ddot{\mathbf{u}}_{t}) \bigr)\\ &\qquad- C_4 \ddot{\mathbf{u}}_{t} - C_5 \dot{\mathbf{u}}_{t} - C_6 \mathbf{u}_{t}\\ \end{aligned}

A\mathbf{A}ut\mathbf{u}_{t} の組み立て

A\mathbf{A}ut\mathbf{u}_{t} の組み立て方について2次元の例を作って見ていきます。 4×44 \times 4uu を考えます。

Image of gibbs phenomenon.

そのまま2次元配列として実装できますが、ソルバで解くときだけ1次元配列にする必要があります。次の図では横書きの文章を書くように、左から右に進み、右端に辿りついたら次の行の左端からまた進むというようにインデックスを振っています。

Image of gibbs phenomenon.
for (var x = 0; x < xLength; ++x) { for (var y = 0; y < yLength; ++y) { u_array[x + y * xLength] = u[x][y] } }

これで ut\mathbf{u}_{t} が組み立てられました。

A\mathbf{A} を組み立てます。整理した式の左辺にある関数 neighbor\mathtt{neighbor} を展開して2次元の形にします。このときインデックスが配列の端から外に出ないよう例外処理が必要です。

C0u¨t+1+C1neighbor(u¨t+1)=C0u¨t+1+C1iaxis(n)(u¨t,i1+u¨t,i+1)=C0u¨t+1+C1(u¨t,x1+u¨t,x+1)+C1(u¨t,y1+u¨t,y+1)=C0u¨t+1+C1sum_x(u¨t+1)+C1sum_y(u¨t+1) \begin{aligned} C_0 \ddot{u}_{t+1} + C_1 \,\mathtt{neighbor}(\ddot{u}_{t+1}) &= C_0 \ddot{u}_{t+1} + C_1 \sum_{i\,\in \mathtt{axis}(n)} ( \ddot{u}_{t,i-1} + \ddot{u}_{t,i+1})\\ &= C_0 \ddot{u}_{t+1} + C_1 (\ddot{u}_{t,x-1} + \ddot{u}_{t,x+1}) + C_1 (\ddot{u}_{t,y-1} + \ddot{u}_{t,y+1}) \\ &= C_0 \ddot{u}_{t+1} + C_1 \mathtt{sum\_x}(\ddot{u}_{t+1}) + C_1 \mathtt{sum\_y}(\ddot{u}_{t+1})\\ \end{aligned}
sum_x(ut)={but,x1,x=0but,x+1,x=nx1ut,x1+ut,x+1,otherwisesum_y(ut)={but,y1,y=0but,y+1,y=ny1ut,y1+ut,y+1,otherwise \begin{aligned} \mathtt{sum\_x}(u_{t}) &= \begin{cases} b u_{t,x - 1} &,\,x=0\\ b u_{t,x + 1} &,\,x=n_x - 1\\ u_{t,x-1} + u_{t,x+1} &,\,\text{otherwise} \end{cases}\\ \mathtt{sum\_y}(u_{t}) &= \begin{cases} b u_{t,y-1} &,\,y=0\\ b u_{t,y+1} &,\,y=n_y - 1\\ u_{t,y-1} + u_{t,y+1} &,\,\text{otherwise} \end{cases}\\ \end{aligned}

nxn_xnyn_y はそれぞれ xx 軸と yy 軸のインデックスの数です。

境界条件は bb の値で変更できます。 b=1b=1 のとき固定端、 b=2b=2 のとき自由端になります。

4×44 \times 4 のときの A\mathbf{A} を組み立てます。

A=[C0bC100bC100000000000C1C0C100bC100000000000C1C0C100bC100000000000bC1C0000bC100000000C1000C0bC100C100000000C100C1C0C100C100000000C100C1C0C100C100000000C100bC1C0000C100000000C1000C0bC100C100000000C100C1C0C100C100000000C100C1C0C100C100000000C100bC1C0000C100000000bC1000C0bC100000000000bC100C1C0C100000000000bC100C1C0C100000000000bC100bC1C0] \mathbf{A} = \begin{bmatrix} C_0 & b C_1 & 0 & 0 & b C_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ C_1 & C_0 & C_1 & 0 & 0 & b C_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & C_1 & C_0 & C_1 & 0 & 0 & b C_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & b C_1 & C_0 & 0 & 0 & 0 & b C_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ C_1 & 0 & 0 & 0 & C_0 & b C_1 & 0 & 0 & C_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & C_1 & 0 & 0 & C_1 & C_0 & C_1 & 0 & 0 & C_1 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & C_1 & 0 & 0 & C_1 & C_0 & C_1 & 0 & 0 & C_1 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & C_1 & 0 & 0 & b C_1 & C_0 & 0 & 0 & 0 & C_1 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & C_1 & 0 & 0 & 0 & C_0 & b C_1 & 0 & 0 & C_1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & C_1 & 0 & 0 & C_1 & C_0 & C_1 & 0 & 0 & C_1 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & C_1 & 0 & 0 & C_1 & C_0 & C_1 & 0 & 0 & C_1 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & C_1 & 0 & 0 & b C_1 & C_0 & 0 & 0 & 0 & C_1\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & b C_1 & 0 & 0 & 0 & C_0 & b C_1 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & b C_1 & 0 & 0 & C_1 & C_0 & C_1 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & b C_1 & 0 & 0 & C_1 & C_0 & C_1\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & b C_1 & 0 & 0 & b C_1 & C_0\\ \end{bmatrix}

何となく規則性が見えます。

B=[C0bC100C1C0C100C1C0C100bC1C0],D=[C10000C10000C10000C1],0=[0000000000000000] \mathbf{B} = \begin{bmatrix} C_0 & b C_1 & 0 & 0\\C_1 & C_0 & C_1 & 0\\0 & C_1 & C_0 & C_1\\0 & 0 & b C_1 & C_0\\ \end{bmatrix} ,\quad \mathbf{D} = \begin{bmatrix} C_1 & 0 & 0 & 0\\0 & C_1 & 0 & 0\\0 & 0 & C_1 & 0\\0 & 0 & 0 & C_1\\ \end{bmatrix} ,\quad \mathbf{0} = \begin{bmatrix} 0 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\\ \end{bmatrix}
A=[BbD00DBD00DBD00bDB] \mathbf{A} = \begin{bmatrix} \mathbf{B} & b\mathbf{D} & \mathbf{0} & \mathbf{0}\\ \mathbf{D} & \mathbf{B} & \mathbf{D} & \mathbf{0}\\ \mathbf{0} & \mathbf{D} & \mathbf{B} & \mathbf{D}\\ \mathbf{0} & \mathbf{0} & b\mathbf{D} & \mathbf{B}\\ \end{bmatrix}

格子の大きさが nx×nyn_x \times n_y のとき、 B,D,0\mathbf{B},\,\mathbf{D},\,\mathbf{0}nx×nxn_x \times n_x の正方行列、 A\mathbf{A}ny×nyn_y \times n_y の正方行列になります。

デモと実装

ソルバにガウス-ザイデル法を使っています。計算が重たいので反復回数の上限を8にしています。

キャンバスをクリックすると波が起きます。utu_{t} の値が大きいと明るく、小さいと暗く表示されます。赤は ut>1.0u_{t} > 1.0 、緑は ut<1.0u_{t} < -1.0 を表しています。

実装はややこしくなっています。JavaScriptでは参照が直接書けないので array = [{value: value}, ...] というように Object に値を格納することで配列の操作を楽にしています。

実装を読む (github.io)

その他

ガウス-ザイデル法による歪みとインデックスの振り方

ガウス-ザイデル法の反復回数が少ないときに波が歪みます。この歪みは1次元配列への変換時に蛇行するようにインデックスを振ることで改善できますが、端のノードについての例外処理が増えます。

Image of gibbs phenomenon.
for (var x = 0; x < xLength; ++x) { for (var y = 0; y < yLength; ++y) { var index = y % 2 == 0 ? x + y * xLength : (y + 1) * xLength - x - 1 u_array[index] = u[x][y] } }

ソルバの発散を防ぐ

2回の反復を1セットにして、1回目はインデックスを昇順、2回目は降順にたどって計算することで、反復回数が少ないときでもソルバによる発散を防げます。

utu_{t}5×35 \times 3 のときの A\mathbf{A}

ut\mathbf{u}_{t} にインデックスを振ります。

Image of gibbs phenomenon.

A\mathbf{A} を組み立てます。

A=[C0bC1000bC1000000000C1C0C1000bC1000000000C1C0C1000bC1000000000C1C0C1000bC1000000000bC1C00000bC100000C10000C0bC1000C100000C1000C1C0C1000C100000C1000C1C0C1000C100000C1000C1C0C1000C100000C1000bC1C00000C100000bC10000C0bC1000000000bC1000C1C0C1000000000bC1000C1C0C1000000000bC1000C1C0C1000000000bC1000bC1C0] \mathbf{A} = \begin{bmatrix} C_0 & b C_1 & 0 & 0 & 0 & b C_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ C_1 & C_0 & C_1 & 0 & 0 & 0 & b C_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & C_1 & C_0 & C_1 & 0 & 0 & 0 & b C_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & C_1 & C_0 & C_1 & 0 & 0 & 0 & b C_1 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & b C_1 & C_0 & 0 & 0 & 0 & 0 & b C_1 & 0 & 0 & 0 & 0 & 0\\ C_1 & 0 & 0 & 0 & 0 & C_0 & b C_1 & 0 & 0 & 0 & C_1 & 0 & 0 & 0 & 0\\ 0 & C_1 & 0 & 0 & 0 & C_1 & C_0 & C_1 & 0 & 0 & 0 & C_1 & 0 & 0 & 0\\ 0 & 0 & C_1 & 0 & 0 & 0 & C_1 & C_0 & C_1 & 0 & 0 & 0 & C_1 & 0 & 0\\ 0 & 0 & 0 & C_1 & 0 & 0 & 0 & C_1 & C_0 & C_1 & 0 & 0 & 0 & C_1 & 0\\ 0 & 0 & 0 & 0 & C_1 & 0 & 0 & 0 & b C_1 & C_0 & 0 & 0 & 0 & 0 & C_1\\ 0 & 0 & 0 & 0 & 0 & b C_1 & 0 & 0 & 0 & 0 & C_0 & b C_1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & b C_1 & 0 & 0 & 0 & C_1 & C_0 & C_1 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & b C_1 & 0 & 0 & 0 & C_1 & C_0 & C_1 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & b C_1 & 0 & 0 & 0 & C_1 & C_0 & C_1\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & b C_1 & 0 & 0 & 0 & b C_1 & C_0\\ \end{bmatrix}

まとめます。

B=[C0bC1000C1C0C1000C1C0C1000C1C0C1000bC1C0],D=[C100000C100000C100000C100000C1],0=[0000000000000000000000000] \mathbf{B} = \begin{bmatrix} C_0 & b C_1 & 0 & 0 & 0\\ C_1 & C_0 & C_1 & 0 & 0\\ 0 & C_1 & C_0 & C_1 & 0\\ 0 & 0 & C_1 & C_0 & C_1\\ 0 & 0 & 0 & b C_1 & C_0\\ \end{bmatrix} ,\quad \mathbf{D} = \begin{bmatrix} C_1 & 0 & 0 & 0 & 0\\ 0 & C_1 & 0 & 0 & 0\\ 0 & 0 & C_1 & 0 & 0\\ 0 & 0 & 0 & C_1 & 0\\ 0 & 0 & 0 & 0 & C_1\\ \end{bmatrix} ,\quad \mathbf{0} = \begin{bmatrix} 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0\\ \end{bmatrix}
A=[BbD0DBD0bDB] \mathbf{A} = \begin{bmatrix} \mathbf{B} & b\mathbf{D} & \mathbf{0}\\ \mathbf{D} & \mathbf{B} & \mathbf{D}\\ \mathbf{0} & b\mathbf{D} & \mathbf{B}\\ \end{bmatrix}

nn 次元の A\mathbf{A}

横書き順インデックスを使うときの nn 次元の A\mathbf{A} を組み立てます。 4×4×4×4 \times 4 \times 4 \times \dots の場合を考えてみます。

B=[C0bC100C1C0C100C1C0C100bC1C0]B2=[BbD00DBD00DBD00bDB]B3=[B2bD20202D2B2D20202D2B2D20202bD2B2]D=[C10000C10000C10000C1]D2=[D0000D0000D0000D]D3=[D202020202D202020202D202020202D2]0=[0000000000000000]02=[0000000000000000]03=[02020202020202020202020202020202] \begin{aligned} \begin{aligned} \mathbf{B} &= \begin{bmatrix} C_0 & b C_1 & 0 & 0\\ C_1 & C_0 & C_1 & 0\\ 0 & C_1 & C_0 & C_1\\ 0 & 0 & b C_1 & C_0\\ \end{bmatrix} \\ \mathbf{B}^{2} &= \begin{bmatrix} \mathbf{B} & b \mathbf{D} & \mathbf{0} & \mathbf{0}\\ \mathbf{D} & \mathbf{B} & \mathbf{D} & \mathbf{0}\\ \mathbf{0} & \mathbf{D} & \mathbf{B} & \mathbf{D}\\ \mathbf{0} & \mathbf{0} & b \mathbf{D} & \mathbf{B}\\ \end{bmatrix} \\ \mathbf{B}^{3} &= \begin{bmatrix} \mathbf{B}^{2} & b \mathbf{D}^{2} & \mathbf{0}^{2} & \mathbf{0}^{2}\\ \mathbf{D}^{2} & \mathbf{B}^{2} & \mathbf{D}^{2} & \mathbf{0}^{2}\\ \mathbf{0}^{2} & \mathbf{D}^{2} & \mathbf{B}^{2} & \mathbf{D}^{2}\\ \mathbf{0}^{2} & \mathbf{0}^{2} & b \mathbf{D}^{2} & \mathbf{B}^{2}\\ \end{bmatrix} \\ &\hspace{6.5em}\vdots \end{aligned} \quad \begin{aligned} \mathbf{D} &= \begin{bmatrix} C_1 & 0 & 0 & 0\\ 0 & C_1 & 0 & 0\\ 0 & 0 & C_1 & 0\\ 0 & 0 & 0 & C_1\\ \end{bmatrix} \\ \mathbf{D}^{2} &= \begin{bmatrix} \mathbf{D} & \mathbf{0} & \mathbf{0} & \mathbf{0}\\ \mathbf{0} & \mathbf{D} & \mathbf{0} & \mathbf{0}\\ \mathbf{0} & \mathbf{0} & \mathbf{D} & \mathbf{0}\\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{D}\\ \end{bmatrix} \\ \mathbf{D}^{3} &= \begin{bmatrix} \mathbf{D}^{2} & \mathbf{0}^{2} & \mathbf{0}^{2} & \mathbf{0}^{2}\\ \mathbf{0}^{2} & \mathbf{D}^{2} & \mathbf{0}^{2} & \mathbf{0}^{2}\\ \mathbf{0}^{2} & \mathbf{0}^{2} & \mathbf{D}^{2} & \mathbf{0}^{2}\\ \mathbf{0}^{2} & \mathbf{0}^{2} & \mathbf{0}^{2} & \mathbf{D}^{2}\\ \end{bmatrix} \\ &\hspace{6em}\vdots \end{aligned} \quad \begin{aligned} \mathbf{0} &= \begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ \end{bmatrix} \\ \mathbf{0}^{2} &= \begin{bmatrix} \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0}\\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0}\\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0}\\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0}\\ \end{bmatrix} \\ \mathbf{0}^{3} &= \begin{bmatrix} \mathbf{0}^{2} & \mathbf{0}^{2} & \mathbf{0}^{2} & \mathbf{0}^{2}\\ \mathbf{0}^{2} & \mathbf{0}^{2} & \mathbf{0}^{2} & \mathbf{0}^{2}\\ \mathbf{0}^{2} & \mathbf{0}^{2} & \mathbf{0}^{2} & \mathbf{0}^{2}\\ \mathbf{0}^{2} & \mathbf{0}^{2} & \mathbf{0}^{2} & \mathbf{0}^{2}\\ \end{bmatrix} \\ &\hspace{5.5em}\vdots \end{aligned} \end{aligned}
A=[BnbDn0n0nDnBnDn0n0nDnBnDn0n0nbDnBn] \mathbf{A} = \begin{bmatrix} \mathbf{B}^{n} & b\mathbf{D}^{n} & \mathbf{0}^{n} & \mathbf{0}^{n}\\ \mathbf{D}^{n} & \mathbf{B}^{n} & \mathbf{D}^{n} & \mathbf{0}^{n}\\ \mathbf{0}^{n} & \mathbf{D}^{n} & \mathbf{B}^{n} & \mathbf{D}^{n}\\ \mathbf{0}^{n} & \mathbf{0}^{n} & b\mathbf{D}^{n} & \mathbf{B}^{n}\\ \end{bmatrix}