📂Numerical Analysis

Parasitic Solution

Definition¹

A parasitic solution, as the name literally suggests a solution that ‘parasitizes’, refers to a term that grows in magnitude and alternates in sign as the method progresses. It helps to picture a sequence like $a_{n} = 2^{-n} + (-2)^{n}$ failing to converge because of $ (-2)^{n}$. Calling such terms ‘parasitic’ is quite an intuitive and apt naming, in that they hinder convergence.

Example: Dahlquist Problem

As an example, consider $\begin{cases} y ' = \lambda y \\ y(0) = 1 \end{cases}$, whose solution is obtained exactly as $Y = e^{ \lambda x}$. However, if what we need are concrete values, we have no choice but to consider numerical methods. Let us use the midpoint method for the computation.

Midpoint method: For a continuous function $f$ defined on $D \subset \mathbb{R}^2$, suppose the initial value problem $\begin{cases} y ' = f(x,y) \\ ( y( x_{0} ) , y (x_{1} )) = (Y_{0} ,Y_{1} ) \end{cases}$ is given. Suppose the interval $(a,b)$ is split into node points such as $a \le x_{0} < x_{1} < \cdots < x_{n} < \cdots x_{N} \le b$. In particular, if we set $x_{j} = x_{0} + j h$ for sufficiently small $h > 0$, then for the initial value $y_{0} = Y_{0}$, $$ y_{n+1} = y_{n-1} + 2 h f ( x_{n} , y_{n} ) $$

Applying the midpoint method to the problem gives the following. $$ y_{n+1} = y_{n-1} + 2h \lambda y_{n} $$ Let us approach this with the method for solving second-order linear homogeneous differential equations, which finds the solution through a specific equation. Assuming $y_{n} = r^{n}$, $$ r^{n+1} = r^{n-1} + 2 h \lambda r^{n} $$ Cancelling $r^{n-1}$ from both sides and rearranging into a quadratic equation, $$ r^2 - 2h \lambda r - 1 = 0 $$ Solving via the quadratic formula, $$ r_{0} = h \lambda + \sqrt{ 1 + h^2 \lambda^2 } $$

$$ r_{1} = h \lambda - \sqrt{ 1 + h^2 \lambda^2 } $$ The general solution, for some $\beta_{0} , \beta_{1}$, is $$ y_{n} = \beta_{0} r_{0}^{n} + \beta_{1} r_{1}^{n} $$ Substituting $n=0,1$, $$ \begin{cases} y_{0} = \beta_{0} + \beta_{1} \\ y_{1} = \beta_{0} r_{0} + \beta_{1} r_{1} \end{cases} $$ Meanwhile, since we already know the exact solution $Y = e^{ \lambda x}$, we can obtain $$ \begin{cases} y_{0} = 1 = \beta_{0} + \beta_{1} \\ y_{1} = e^{ \lambda h } = \beta_{0} r_{0} + \beta_{1} r_{1} \end{cases} $$ Solving this for $\beta_{0}$ and $\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} $$ Expanding $e^{\lambda h}$ with the Maclaurin expansion gives $\displaystyle e^{\lambda h} = 1 + \lambda h + {{\lambda^2 h^2} \over {2}} + O (h^3 \lambda^3)$, so $$ \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*} $$ Expanding $\sqrt{ 1+ h^2 \lambda^2 }$ with the Maclaurin expansion gives $\displaystyle \sqrt{ 1+ h^2 \lambda^2 } = 1 + {{\lambda^2 h^2} \over {2}} + O (h^4 \lambda^4)$, so $$ \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*} $$ Likewise for $\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*} $$

How to move big O notation from the denominator up to the numerator: For $a \ne 0$ and $p>0$, $n \in \mathbb{N}$, $\displaystyle {{1} \over { \sqrt[p]{a + O ( h^n ) } }} = {{1} \over { \sqrt[p]{a } }}+ O(h^n)$

Putting it all together, we obtain $$ \begin{align*} \beta_{0} =& 1+ O (h^{3} \lambda^{3} ) \\ \beta_{1} =& O (h^{3} \lambda^{3} ) \end{align*} $$ That is, as $h \to 0$ we have $\beta_{0} \to 1$ and $\beta_{1} \to 0$, so $$ y_{n} = \beta_{0} r_{0}^{n} + \beta_{1} r_{1}^{n} \to r_{0}^{n} $$ Now the question is what happens to this general solution as $n$ grows while $h>0$ is fixed. If $\lambda > 0$, there is nothing to worry about: since $r_{0} > | r_{1} | > 0$, $\beta_{0} r_{0}^{b}$ grows much faster than $\beta_{1} r_{1}^{n}$. But if $\lambda <0$ the story is different: if $0 < r_{0} < 1$ and $r_{1} < -1$, then $\beta_{1} r_{1}^{n}$ flips its sign every time $n$ increases, while its absolute value comes to overwhelm $\beta_{0} r_{0}^{n}$.

It is precisely this $\beta_{1} r_{1}^{n}$ that we call the parasitic solution, and because of this danger, the midpoint method is said to have weak stability. Therefore, at least when the sign of $\displaystyle { {\partial f(x, Y(x)) } \over {\partial y}}$ is negative, one must check the formulas to see whether such a problem is present.

As an example, looking at the result of solving an initial value problem such as $\begin{cases} y ' = x - y^2 \\ y(0) = 0 \end{cases}$ with the midpoint method, it seems to go well at first, but from the third row on the solution suddenly starts to oscillate.

In almost the same way, 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}) $$ can also be shown to have weak stability. A problem such as $\begin{cases} y ' = \lambda y \\ y(0) = 1 \end{cases}$, favored for this purpose, is called the Dahlquist problem.