📂数値解析

パラサイティック・ソリューション

定義¹

パラサイティック・ソリューション^{parasitic solution}とは、直訳すると「寄生する解」という意味で、メソッドが進むにつれて大きさが増し、符号が入れ替わるなどの項をいう。$a_{n} = 2^{-n} + (-2)^{n}$という数列が$ (-2)^{n}$のせいで収束しないことを想像するとよい。このような項に「パラサイティック」という表現を使うのは、収束を妨げるという点でかなり直感的で良い命名だと言えるだろう。

例: ダールキスト問題

例として$\begin{cases} y ' = \lambda y \\ y(0) = 1 \end{cases}$を考えてみると、その解は$Y = e^{ \lambda x}$と正確に求められる。しかし、我々が必要とするのが具体的な値であれば、数値解析的な方法を考慮せざるを得ない。計算のためにミッドポイントメソッドを使ってみよう。

ミッドポイントメソッド: $D \subset \mathbb{R}^2$で定義された連続関数$f$に対して、初期値問題$\begin{cases} y ' = f(x,y) \\ ( y( x_{0} ) , y (x_{1} )) = (Y_{0} ,Y_{1} ) \end{cases}$が与えられている。区間$(a,b)$を$a \le x_{0} < x_{1} < \cdots < x_{n} < \cdots x_{N} \le b$のようなノードポイントに分割したとしよう。特に十分小さい$h > 0$に対して$x_{j} = x_{0} + j h$とすると、初期値$y_{0} = Y_{0}$に対して $$ y_{n+1} = y_{n-1} + 2 h f ( x_{n} , y_{n} ) $$

ミッドポイントメソッドを問題に適用すると次のようになる。 $$ y_{n+1} = y_{n-1} + 2h \lambda y_{n} $$ 特定の方程式を通じて解を求める2階線形同次微分方程式の解法でアプローチしてみよう。$y_{n} = r^{n}$と仮定してみると $$ r^{n+1} = r^{n-1} + 2 h \lambda r^{n} $$ 両辺から$r^{n-1}$を消去し、2次方程式として整理すると $$ r^2 - 2h \lambda r - 1 = 0 $$ 解の公式で解くと $$ r_{0} = h \lambda + \sqrt{ 1 + h^2 \lambda^2 } $$

$$ r_{1} = h \lambda - \sqrt{ 1 + h^2 \lambda^2 } $$ 一般解はある$\beta_{0} , \beta_{1}$に対して $$ y_{n} = \beta_{0} r_{0}^{n} + \beta_{1} r_{1}^{n} $$ $n=0,1$を代入すると $$ \begin{cases} y_{0} = \beta_{0} + \beta_{1} \\ y_{1} = \beta_{0} r_{0} + \beta_{1} r_{1} \end{cases} $$ 一方、我々はすでに正確な解として$Y = e^{ \lambda x}$を知っているので $$ \begin{cases} y_{0} = 1 = \beta_{0} + \beta_{1} \\ y_{1} = e^{ \lambda h } = \beta_{0} r_{0} + \beta_{1} r_{1} \end{cases} $$ を求めることができる。これを$\beta_{0}$と$\beta_{1}$について解くと $$ \begin{cases} \displaystyle \beta_{0} = {{e^{ \lambda h} - r_{1} } \over {2 \sqrt{ 1+ h^2 \lambda^2 } }} \\ \displaystyle \beta_{1} = {{r_{0} - e^{ \lambda h} } \over {2 \sqrt{ 1+ h^2 \lambda^2 } } } \end{cases} $$ $e^{\lambda h}$をマクローリン展開すると$\displaystyle e^{\lambda h} = 1 + \lambda h + {{\lambda^2 h^2} \over {2}} + O (h^3 \lambda^3)$なので $$ \begin{align*} \displaystyle \beta_{0} =& {{e^{ \lambda h} - r_{1} } \over {2 \sqrt{ 1+ h^2 \lambda^2 } }} \\ =& {{1 + \lambda h + {{\lambda^2 h^2} \over {2}} + O (h^3 \lambda^3 ) - h \lambda + \sqrt{ 1 + h^2 \lambda^2 } } \over {2 \sqrt{ 1+ h^2 \lambda^2 } }} \\ =& {{1 + {{\lambda^2 h^2} \over {2}} - \sqrt{ 1 + h^2 \lambda^2 } + 2 \sqrt{ 1 + h^2 \lambda^2 } + O (h^3 \lambda^3 ) } \over {2 \sqrt{ 1+ h^2 \lambda^2 } }} \\ =& 1 + {{1 + {{\lambda^2 h^2} \over {2}} - \sqrt{ 1 + h^2 \lambda^2 } + O (h^3 \lambda^3 ) } \over {2 \sqrt{ 1+ h^2 \lambda^2 } }} \end{align*} $$ $\sqrt{ 1+ h^2 \lambda^2 }$をマクローリン展開すると$\displaystyle \sqrt{ 1+ h^2 \lambda^2 } = 1 + {{\lambda^2 h^2} \over {2}} + O (h^4 \lambda^4)$なので $$ \begin{align*} \displaystyle \beta_{0} =& 1 + {{1 + {{\lambda^2 h^2} \over {2}} - \sqrt{ 1 + h^2 \lambda^2 } +O (h^3 \lambda^3 ) } \over {2 \sqrt{ 1+ h^2 \lambda^2 } }} \\ =& 1 + {{1 + {{\lambda^2 h^2} \over {2}} + O (h^3 \lambda^3 ) - 1 - {{\lambda^2 h^2} \over {2}} + O (h^4 \lambda^4) } \over {2 \sqrt{ 1+ h^2 \lambda^2 } }} \end{align*} $$ $\beta_{1}$についても同様に $$ \begin{align*} \beta_{1} =& {{ r_{0} - e^{ \lambda h} } \over {2 \sqrt{ 1+ h^2 \lambda^2 } }} \\ =& { { h \lambda + \sqrt{ 1+ h^2 \lambda^2 } - e^{h \lambda } } \over {2 \sqrt{ 1+ h^2 \lambda^2 } }} \\ =& {{ h \lambda + 1 + {{\lambda^2 h^2} \over {2}} + O (h^4 \lambda^4 ) - 1 - h \lambda - {{\lambda^2 h^2} \over {2}} - O (h^3 \lambda^3) } \over {2 \sqrt{ 1+ h^2 \lambda^2 } }} \end{align*} $$

ビッグオー記法が分母にあるとき分子に上げる方法: $a \ne 0$と$p>0$、$n \in \mathbb{N}$に対して$\displaystyle {{1} \over { \sqrt[p]{a + O ( h^n ) } }} = {{1} \over { \sqrt[p]{a } }}+ O(h^n)$

整理すると $$ \begin{align*} \beta_{0} =& 1+ O (h^{3} \lambda^{3} ) \\ \beta_{1} =& O (h^{3} \lambda^{3} ) \end{align*} $$ を得る。すなわち$h \to 0$のとき$\beta_{0} \to 1$であり$\beta_{1} \to 0$なので $$ y_{n} = \beta_{0} r_{0}^{n} + \beta_{1} r_{1}^{n} \to r_{0}^{n} $$ さて問題は、$h>0$が定められているとき、$n$が大きくなるにつれてこの一般解がどうなるかだ。もし$\lambda > 0$なら悩むまでもなく$r_{0} > | r_{1} | > 0$なので、$\beta_{0} r_{0}^{b}$が$\beta_{1} r_{1}^{n}$よりはるかに速く大きくなる。しかし$\lambda <0$なら話は変わり、もし$0 < r_{0} < 1$かつ$r_{1} < -1$であれば、$\beta_{1} r_{1}^{n}$は$n$が増えるたびに符号を変えながら、その絶対値は$\beta_{0} r_{0}^{n}$を圧倒するようになる。

このとき、まさにこの$\beta_{1} r_{1}^{n}$をパラサイティック・ソリューションと呼び、このような危険のためにミッドポイントメソッドは弱い安定性^{weak Stability}を持つという。したがって、少なくとも$\displaystyle { {\partial f(x, Y(x)) } \over {\partial y}}$の符号が負のときは、このような問題がないか数式による確認が必ず必要である。

例として$\begin{cases} y ' = x - y^2 \\ y(0) = 0 \end{cases}$のような初期値問題をミッドポイントメソッドで解いた結果を見ると、最初はうまくいくように見えるが、3行目からは解が突然揺れ動き始めることが確認できる。

これとほとんど同じ方法で、ミルンメソッド^{milne's method}

$$ y_{n+1} = y_{n-1} + {{h} \over {3}} [ f(x_{n-1} , y_{n-1}) + f(x_{n} , y_{n})] + f(x_{n+1} , y_{n+1}) $$ も弱い安定性を持つことが示せる。このとき好んで使われる$\begin{cases} y ' = \lambda y \\ y(0) = 1 \end{cases}$のような問題をダールキスト問題^{dahlquist problem}という。