Section 1.9 Stiff ODEs
. The explicit Euler method and, in general, methods that come from generalization of Newton-Cotes methods, don`t do well with a type of ODEs called stiff. We show below an elementary example of this behavior. Consider the ODE
\begin{equation*}
\dot x = -kx,\;\;x(0)=1
\end{equation*}
for some \(k>0\) and apply the explicit Euler method to it:
\begin{equation*}
x_{n+1}=x_n-hkx_n = (1-hk)x_n.
\end{equation*}
In this case we can find a closed formula for the \(n\)-th iteration:
\begin{equation*}
x_n = (1-hk)^{n},\;\;n=1,2,\dots
\end{equation*}
If \(k\) is large, the sequence \(x_n\text{,}\) unless \(h\) is chosen small enough, diverges and oscillates (its sign will switch at every step). This is exactly the opposite of the qualitative behavior of the exact solution of the ODE:
\begin{equation*}
x(t) = e^{-kt}.
\end{equation*}
This solution converges monotonically to 0 for all \(k>0\text{.}\) This is bad: a "good" method should, at least, follow the same qualitative behavior of the exact solution. The code below, solving the ODE \(\dot x = -15x\text{,}\) shows had badly the algorithm can behave. The situation does not improve much using Heun's method (try it!).
Stability regions. Consider now the general linear IVP
\begin{equation*}
\dot x = Kx,\;\;x(t_0)=x_0,\;\;K\in\mathbb C.
\end{equation*}
We pose the following question: given a numerical method to solve ODEs with a given step size \(h\) and a \(K\) with negative real part, for which values of \(Kh\) does its numerical solution \(x_n\) behave qualitatively similarly to its exact solution \(x(t)=x_0e^{K(t-t_0)}\text{?}\) The set of all \(z=Kh\) in the complex plane for which the numerical solution behaves qualitatively like the exact solution is called the stability region of the method.
Stability region of the explicit Euler method. The argument given at the beginning of the section shows that
\begin{equation*}
x_n = x_0(1+Kh)^n.
\end{equation*}
Hence, we get the same qualitative behavior as long as \(x_n\) remains bounded for all \(n\text{.}\) In turn, this happens if and only if
\begin{equation*}
|1+Kh|\leq1.
\end{equation*}
We plot this region of the complex plane below.
Stability region of the Heun method. Using the same kind of argument above, we find that
\begin{equation*}
x_n = x_0\left(1+Kh+\frac{K^2h^2}{2}\right)^n.
\end{equation*}
Hence, we get the same qualitative behavior as long as \(x_n\) remains bounded for all \(n\text{.}\) In turn, this happens if and only if
\begin{equation*}
\left|1+Kh+\frac{K^2h^2}{2}\right|\leq1.
\end{equation*}
We plot this region of the complex plane below.
Stability region of the RK4 method. Using the same kind of argument above, we find that
\begin{equation*}
x_n = x_0\left(1+Kh+\frac{K^2h^2}{2}+\frac{K^3h^3}{6}+\frac{K^4h^4}{24}\right)^n.
\end{equation*}
Hence, we get the same qualitative behavior as long as \(x_n\) remains bounded for all \(n\text{.}\) In turn, this happens if and only if
\begin{equation*}
\left|1+Kh+\frac{K^2h^2}{2}+\frac{K^3h^3}{6}+\frac{K^4h^4}{24}\right|\leq1.
\end{equation*}
We plot this region of the complex plane below.