Subsections

• 6.1 固有値を求める
• 6.2 「直観的」説明

6 線形安定性解析

さて、まず、なぜ振動が起きるかということだが、これは原理的には差分化さ れた方程式の固有値を求めればわかる。

つまり、 ${\boldmath {u}}_{j}$ を時刻 $t_{j} = j\Delta t$ での数値解をベク トルとして書いたものだとすると、差分近似は、

$${\boldmath {u}}_{j+1} = A{\boldmath {u}}_{j}$$ (8)

という形に書ける。ここで $A$ は$(n-1)$次元正方行列で、その各要素は差分近似の形から決 まる。この場合には、対角要素とその両側一列以外の全ての要素が 0 である3 重対角行列といわれる形をしている。要素の値は

$$a_{ii} = 1-2 \delta\quad (0<i<n)$$ (9)

$$a_{i,i-1} = a_{i-1,i} = \delta$$ (10)

ここで、 $\delta = \Delta t/\Delta x^2$ である。

この行列は実対称行列なので固有値分解ができて、適当なユニタリ行列$U$を持っ てきて

$$A = U^{-1}\Lambda U$$ (11)

と書ける。ここで$\Lambda$ は固有値行列、つまり対角要素が固有値でそれ以 外が 0の行列である。

このように分解できるので、 ${\boldmath {v}}= U{\boldmath {u}}$ とすれば最初の差分近似式は

$${\boldmath {v}}_{j+1} = \Lambda {\boldmath {v}}_j$$ (12)

つまり、

$$v_{i,j+1} = \lambda_i {\boldmath {v}}_{ij}$$ (13)

という $n-1$ 個の独立な方程式になって、解も自明に求まる。いうまでもな いが、 $U$ の各列は A の固有ベクトルになっている。つまり、$U$の $i$ 行 目を ${\boldmath {x}}_i$ と書くと、

$$\lambda {\boldmath {x}}_i = A{\boldmath {x}}_i$$ (14)

である。

不安定性が起きないための条件は、全ての固有値 $\lambda_i$ の絶対値が $1$ を超えないということである。

6.1 固有値を求める

さて、そういうわけで、固有値$\lambda$がどうなっているか調べてみよう。 式 (14) をもともとの差分近似に入れてちょっと変形する と、

$$(\lambda + 2 \delta -1)u_i = \delta(u_{i+1}+u_{i-1})$$ (15)

となる。 $x$ がいつのまにか $u$ に戻っているが深い意味はない。ここで $\delta = \Delta t/\Delta x^2$ である。さらに、

$$\alpha = \frac{\lambda + 2 \delta -1}{\delta}$$ (16)

とおいてもうちょっと変形すると、結局

$$u_{i+1} - \alpha u_i + u_{i-1} = 0$$ (17)

となる。

断るまでもないと思うが、これは元の偏微分を変数分離して空間方向の関数に ついての常微分方程式を導くのとほとんど同じ操作になっている。これを $u_i$ についての差分方程式と見て解を求めることを考える。

線形差分方程式なので、解は $u_i = Cp^i$ の形である。これを代入すると

$$p^2 - \alpha p + 1 = 0$$ (18)

解は

$$p=\frac{\alpha \pm \sqrt{\alpha^2 -4}}{2}$$ (19)

となる。

真面目にやるには境界条件を満たす固有ベクトルをもとめないといけないが、 面倒なので無限遠境界の場合を考える。この時、 $p$ が実数のものはどち らかで無限大に発散するのでよろしくないので、複素数の場合を考える。この 時は $\vert p\vert = 1$ になっているので無限遠でも発散しない。

なお、有限の固定境界の場合も結局境界条件を満たすような解を作るためには $p$ が複素数でないといけないことがすぐにわかる。このため、 $\vert \alpha \vert < 2$ だけを考えればよい。

このとき、式(16)は

$$\lambda = 1- (2-\alpha) \delta$$ (20)

と書き直せる。$\delta > 0$ なので、$\vert \alpha \vert < 2$ なるある $\alpha$ に ついて$\vert\lambda\vert<1$ であるためには、結局

$$\delta < \frac{2}{2-\alpha}$$ (21)

を満たせば良く、任意の$\alpha$ についてなり立つためには

$$\delta < \frac{1}{2}$$ (22)

つまり

$$\frac{\Delta t}{\Delta x^2} < \frac{1}{2}$$ (23)

であればいいことになる。

上でやったことは、結局空間方向をフーリエ級数展開して、各空間波長に対す る時間発展が安定である（減衰していく）ことを要求しているのと同じである。 このようなやり方を von Neumann の方法による安定性解析という。

6.2 「直観的」説明

上の議論からわかるように、不安定条件に関係する $\alpha$ の値は実際には $\alpha = -2$ だけで、それ以外の $\alpha$ からはもっと緩い条件しかでな い。 $\alpha = -2$ は $p=-1$ に対応するので、これは結局 $u_i$ と $u_{i+1}$ で値が逆転するようなパターンである。そのようなパターン（モー ド）が 1 ステップ計算するとどうなるかをもう一度もとの差分式に戻って書き直してみ ると、結局

$$u_{i,j+1} = (1-4\Delta t /\Delta x^2)u_{ij}$$ (24)

となる。括弧内の絶対値が 1 より大きいと、パターンがどんどん成長していっ てめちゃくちゃな答になるわけである。括弧内の絶対値が 1 より小さいため には、

$$\frac{\Delta t}{\Delta x^2} < \frac{1}{2}$$ (25)

であればよく、前節の結果と同じである。このように、偏微分方程式の場合に は大抵もっとも波長の短いモードが減衰できるかどうかで安定性が決まる。