logo

パラサイティック・ソリューション 📂数値解析

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

定義 1

パラサイティックソリューションparasitic solutionとは、方法が進むにつれてその大きさが大きくなり、符号が変わるような項のことだ。an=2n+(2)na_{n} = 2^{-n} + (-2)^{n}のような数列が(2)n (-2)^{n}のせいで収束しないことを想像するといい。「パラサイティック」という表現は、収束を妨げるという点で直感的で、非常に適切な命名だと言えるだろう。

例:ダールキスト問題

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

ミッドポイントメソッドDR2D \subset \mathbb{R}^2によって定義された連続関数ffに対して、初期値問題{y=f(x,y)(y(x0),y(x1))=(Y0,Y1)\begin{cases} y ' = f(x,y) \\ ( y( x_{0} ) , y (x_{1} )) = (Y_{0} ,Y_{1} ) \end{cases}が与えられたとする。区間(a,b)(a,b)ax0<x1<<xn<xNba \le x_{0} < x_{1} < \cdots < x_{n} < \cdots x_{N} \le bのようなノードポイントで切り分けたとしよう。特に、十分に小さいh>0h > 0に対してxj=x0+jhx_{j} = x_{0} + j hとすると、初期値y0=Y0y_{0} = Y_{0}に対して yn+1=yn1+2hf(xn,yn) y_{n+1} = y_{n-1} + 2 h f ( x_{n} , y_{n} )

ミッドポイントメソッドを問題に適用すると、次のようになる。 yn+1=yn1+2hλyn y_{n+1} = y_{n-1} + 2h \lambda y_{n} 2次線形同次微分方程式の解法を通じて、yn=rny_{n} = r^{n}と仮定しよう。 rn+1=rn1+2hλrn r^{n+1} = r^{n-1} + 2 h \lambda r^{n} 両辺からrn1r^{n-1}を除去し、2次方程式に整理すると r22hλr1=0 r^2 - 2h \lambda r - 1 = 0 解の公式を通じて解くと r0=hλ+1+h2λ2 r_{0} = h \lambda + \sqrt{ 1 + h^2 \lambda^2 }

r1=hλ1+h2λ2 r_{1} = h \lambda - \sqrt{ 1 + h^2 \lambda^2 } 一般解は、あるβ0,β1\beta_{0} , \beta_{1}に対して yn=β0r0n+β1r1n y_{n} = \beta_{0} r_{0}^{n} + \beta_{1} r_{1}^{n} n=0,1n=0,1を代入すると {y0=β0+β1y1=β0r0+β1r1 \begin{cases} y_{0} = \beta_{0} + \beta_{1} \\ y_{1} = \beta_{0} r_{0} + \beta_{1} r_{1} \end{cases} 一方で、すでに正確な解としてY=eλxY = e^{ \lambda x}を知っているため {y0=1=β0+β1y1=eλh=β0r0+β1r1 \begin{cases} y_{0} = 1 = \beta_{0} + \beta_{1} \\ y_{1} = e^{ \lambda h } = \beta_{0} r_{0} + \beta_{1} r_{1} \end{cases} を導き出すことができる。これをβ0\beta_{0}β1\beta_{1}に対して解くと {β0=eλhr121+h2λ2β1=r0eλh21+h2λ2 \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λhe^{\lambda h}マクローリン展開すると、eλh=1+λh+λ2h22+O(h3λ3)\displaystyle e^{\lambda h} = 1 + \lambda h + {{\lambda^2 h^2} \over {2}} + O (h^3 \lambda^3)であるため β0=eλhr121+h2λ2=1+λh+λ2h22+O(h3λ3)hλ+1+h2λ221+h2λ2=1+λ2h221+h2λ2+21+h2λ2+O(h3λ3)21+h2λ2=1+1+λ2h221+h2λ2+O(h3λ3)21+h2λ2 \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*} 同様に1+h2λ2\sqrt{ 1+ h^2 \lambda^2 }をマクローリン展開すると、1+h2λ2=1+λ2h22+O(h4λ4)\displaystyle \sqrt{ 1+ h^2 \lambda^2 } = 1 + {{\lambda^2 h^2} \over {2}} + O (h^4 \lambda^4)であるため β0=1+1+λ2h221+h2λ2+O(h3λ3)21+h2λ2=1+1+λ2h22+O(h3λ3)1λ2h22+O(h4λ4)21+h2λ2 \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*} β1\beta_{1}に対しても同じことが言える。 β1=r0eλh21+h2λ2=hλ+1+h2λ2ehλ21+h2λ2=hλ+1+λ2h22+O(h4λ4)1hλλ2h22O(h3λ3)21+h2λ2 \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*}

分母にビッグオーノーテーションがある場合、分子に移す方法a0a \ne 0p>0p>0nNn \in \mathbb{N}に対して1a+O(hn)p=1ap+O(hn)\displaystyle {{1} \over { \sqrt[p]{a + O ( h^n ) } }} = {{1} \over { \sqrt[p]{a } }}+ O(h^n)を適用する。

結論として、 β0=1+O(h3λ3)β1=O(h3λ3) \begin{align*} \beta_{0} =& 1+ O (h^{3} \lambda^{3} ) \\ \beta_{1} =& O (h^{3} \lambda^{3} ) \end{align*} を得る。つまり、h0h \to 0の時、β01\beta_{0} \to 1であり、β10\beta_{1} \to 0なので yn=β0r0n+β1r1nr0n y_{n} = \beta_{0} r_{0}^{n} + \beta_{1} r_{1}^{n} \to r_{0}^{n} 今、問題は、nnが固定されたときに、λ>0\lambda > 0が大きくなると、この一般解がどうなるかである。もしλ>0\lambda > 0なら、心配することはなく、r0>r1>0r_{0} > | r_{1} | > 0により、β0r0b\beta_{0} r_{0}^{b}β1r1n\beta_{1} r_{1}^{n}よりもずっと早く大きくなる。しかし、λ<0\lambda <0なら、話は変わる。もし0<r0<10 < r_{0} < 1であり、r1<1r_{1} < -1なら、nnが増加するたびにβ1r1n\beta_{1} r_{1}^{n}は符号を変えながら、その絶対値がβ0r0n\beta_{0} r_{0}^{n}を圧倒することになる。

この時、β1r1n\beta_{1} r_{1}^{n}パラサイティックソリューションと呼び、この危険のためにミッドポイントメソッドが弱い安定性weak Stabilityを持つと言われる。したがって、少なくともf(x,Y(x))y\displaystyle { {\partial f(x, Y(x)) } \over {\partial y}}の符号が負である場合、この問題がないか数学的に確認することが絶対に必要である。

20181009\_124502.png

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

これとほぼ同じ方法で、ミルンメソッドmilne’s method

yn+1=yn1+h3[f(xn1,yn1)+f(xn,yn)]+f(xn+1,yn+1) 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}) も弱い安定性を持つことが示される。特に好んで使われる{y=λyy(0)=1\begin{cases} y ' = \lambda y \\ y(0) = 1 \end{cases}と同じ問題をダールキスト問題dahlquist Problemと言う。


  1. Atkinson. (1989). An Introduction to Numerical Analysis(2nd Edition): p364~365. ↩︎