ディリクレ境界条件が与えられた熱方程式の初期値問題に対する数値解析的解法
##例
$$ \begin{cases} u_{t} = \gamma u_{xx} \\ u(t,0) = u(t,l) = 0 \\ u(0,x) = f(x) \end{cases} $$
与えられた問題には代数的解があるほど簡単で、しかし、微分方程式を解くための数値解析法をなぜ学ぶのかを明快に示す例にもなる。これは単に$y ' = f(x,y)$形式の微分方程式を解くことが、偏微分方程式の解法にも繋がるからだ。
解法
代数的解法とは異なり、$\gamma$の値は重要ではないので、便宜上$\gamma = 1$としよう。
$$ x_{j} = \delta j $$ 直線空間line space$x$のステップサイズを$\delta$とすると、$j= 0, 1, \cdots , m$に対して上記のように表せる。
$x$に対する二階微分の近似として $$ U_{xx} \approx {{ U(t, x_{j+1}) - 2 U(t, x_{j}) + U(t, x_{j-1}) } \over { \delta^2 }} $$ これで$t$に関する方程式 $$ u_{t} (t, x_{j}) = {{ u(t, x_{j+1}) - 2 u(t, x_{j}) + u(t, x_{j-1}) } \over { \delta^2 }} $$ を$m-1$個得ることができる。これを行列で表すと $$ \begin{bmatrix} u_{t} (t, x_{1}) \\ u_{t} (t, x_{2}) \\ u_{t} (t, x_{3}) \\ \vdots \\ u_{t} (t, x_{m-2}) \\ u_{t} (t, x_{m-1}) \end{bmatrix} = {{1} \over {\delta^2}} \begin{bmatrix} -2 & 1 & 0 & \cdots & 0 & 0 \\ 1 & -2 & 1 & \cdots & 0 & 0 \\ 0 & 1 & -2 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -2 & 1 \\ 0 & 0 & 0 & \cdots & 1 & -2 \end{bmatrix} \begin{bmatrix} u (t, x_{1}) \\ u (t, x_{2}) \\ u (t, x_{3}) \\ \vdots \\ u (t, x_{m-2}) \\ u (t, x_{m-1}) \end{bmatrix} $$ であり、ページを節約するために $$ \mathbf{u}’(t) :=\begin{bmatrix} u_{t} (t, x_{1}) \\ u_{t} (t, x_{2}) \\ u_{t} (t, x_{3}) \\ \vdots \\ u_{t} (t, x_{m-2}) \\ u_{t} (t, x_{m-1}) \end{bmatrix} $$
$$ \Lambda := {{1} \over {\delta^2}} \begin{bmatrix} -2 & 1 & 0 & \cdots & 0 & 0 \\ 1 & -2 & 1 & \cdots & 0 & 0 \\ 0 & 1 & -2 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -2 & 1 \\ 0 & 0 & 0 & \cdots & 1 & -2 \end{bmatrix} $$
$$ \mathbf{u} (t) := \begin{bmatrix} u (t, x_{1}) \\ u (t, x_{2}) \\ u (t, x_{3}) \\ \vdots \\ u (t, x_{m-2}) \\ u (t, x_{m-1}) \end{bmatrix} $$ とし、$\mathbf{u} ’ (t) = \Lambda \mathbf{u} (t)$となる。今、時間$t$のステップサイズを$h$とし、$t_{k} = hk$として、いくつかの方法を適用して数値解を求めればよい。
■
教訓: スティッフな微分方程式
しかし、これはあくまで理論上の話で、実際に数値解を求めるには、かなりの困難が伴う。
$\Lambda$の固有値を知るために、$\Lambda - \lambda I_{m-1}$の行列式を求めてみよう。
三重対角行列の行列式: 三重対角行列$X_{n} := \begin{bmatrix} x & 1 & 0 & \cdots & 0 & 0 \\ 1 & x & 1 & \cdots & 0 & 0 \\ 0 & 1 & x & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & x & 1 \\ 0 & 0 & 0 & \cdots & 1 & x \end{bmatrix}$について、$\displaystyle | x_{n}| = U_{n} \left( {{x} \over {2}} \right)$
$$ \Lambda - \lambda I_{m-1} = {{1} \over { \delta^2 }} \begin{bmatrix} -2 - \delta^2 \lambda & 1 & 0 & \cdots & 0 & 0 \\ 1 & -2 - \delta^2 & 1 & \cdots & 0 & 0 \\ 0 & 1 & -2 - \delta^2 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -2 - \delta^2 & 1 \\ 0 & 0 & 0 & \cdots & 1 & -2 - \delta^2 \end{bmatrix} $$ 三重対角行列である$(m - 1) \times (m-1)$なので、行列式は $$ \det \left( \Lambda - \lambda I_{m-1} \right) = {{1} \over { \delta^2 }} U_{m-1} \left( {{-2 - \delta^2 \lambda } \over {2}} \right) $$ のように求められる。ここで、$U_{m-1}$は$m-1$次の第二種チェビシェフ多項式を意味する。
- [1]: $\displaystyle U_{n} (X)$の根は$\displaystyle x_{k} = \cos \left( {{k} \over {n+1}} \pi \right)$、$k=1, \cdots , n$
- [2]: $U_{n} (-x) = (-1)^{n} U_{n} (x)$
括弧から$-1$を取り出せば $$ {{(-1)^{m-1}} \over { \delta^2 }} U_{m-1} \left( {{2 + \delta^2 \lambda } \over {2}} \right) $$ であり、固有値を求めるためには $$ U_{m-1} \left( {{2 + \delta^2 \lambda } \over {2}} \right) $$ だけ気にすれば十分である。$\displaystyle U_{m-1} \left( {{2 + \delta^2 \lambda } \over {2}} \right)$の$j$番目の根は$j=1, \cdots , (m-1)$に対して $$ {{2 + \delta^2 \lambda_{j}} \over {2}} = \cos \left( {{j} \over {m}} \pi \right) $$ であり、両辺に$2$を掛けると $$ \begin{align*} \displaystyle 2+\delta^2 \lambda_{j} =& 2 \cos \left( {{j} \over {m}} \pi \right) \\ =& 2 \cos \left( 2 \cdot {{j} \over {2m}} \pi \right) \\ =& 2 \left[ 1 - 2 \sin^2 \left( {{ j \pi } \over { 2m }} \right) \right] \end{align*} $$ したがって、固有値は $$ \lambda_{j} = - {{4} \over { \delta^2 }} \sin^2 \left( {{j \pi } \over {2m}} \right) $$ のように求められる。しかし、 $$ \begin{align*} \left| \lambda_{m-1} \right| =& \left| {{4} \over { \delta^2 }} \sin^2 \left( {{(m-1) \pi } \over {2m}} \right) \right| \\ \approx & {{4} \over {\delta^2}} \\ =& 4 m^2 \end{align*} $$
であるため、$\displaystyle m = {{1} \over {\delta}}$が大きくなるほど、固有値の大きさも莫大になる。数値解をできるだけ正確に知るために$\delta$を小さくして$[0,1]$をたくさん分割した場合、計算は非常に難しくなる。たとえば$m=100$にしても$\mathbf{u} ’ (t) = \Lambda \mathbf{u} (t)$はスティッフstiffな問題となり、数値的に解くのはかなり困難である。