Using RK4 on the van der Pol oscillator

Question

The van der Pol oscillator is modeled by

$$\frac{d^2y}{dt^2} + \mu(y^2 - 1) \frac{dy}{dt} + y = 0$$

This can be written as a system of first order equations, $$ \begin{pmatrix} \dot x\\ \dot y\end{pmatrix} = \begin{pmatrix} \mu(1-y^2)x - y\\ x\end{pmatrix}, $$

where $x := \dot y$. Consider the case where $y(0) = 0$ and $x(0) = a$, where $a$ is some real number. A numerical approximation for the solution can be obtained using the RK4 method. I have two questions:

1. How do you find the region of stability for the RK4 method?

2. How can you verify RK4 is 4th order accurate? Surely you would need to know the true solution in order to verify this?

Answer 1

$\newcommand\rx{\mathfrak{x}}\newcommand\ru{\mathfrak{u}}\newcommand\mA{\mathfrak{A}}$

1. ) For the definition of the most common stability properties one refers to the properties of solutions of linear systems of ODE. Apply the RK4 method to some linear system $\ru'=\mA\ru$ (for instance the linearization of a non-linear ODE system around a fixed point). You should get exactly the first 5 terms (up to degree 4) of the exponential series.

   The stability condition is that this propagator has to be non-expanding if the spectrum of $\mA$ is contained in the negative half-plane of the complex plane. On the level of singular eigenvalues $\lambda$ and $z=λh$ with step size $h$ this means that $z$ is in the stability region if $$|1+z+z^2/2+z^3/6+z^4/24|\le1.$$

   The half-circle with $\Im(z)\le 0$, $|z|\le 2.5$ is contained in that set.

   If $L$ is a Lipschitz constant (or $L=\|\mA\|$), then this gives a bound $h<2.5/L$. At the upper bound the numerical solution is merely not totally wrong, for a useful solution one needs a much smaller step size, $Lh\in[10^{-3},10^{-1}]$. Step sizes much smaller than that lead to an accumulation of floating point errors that dominate the method error.

   See Maximum timestep for RK4 for visualizations of the stability region and the effect of too small or too large step sizes.

2. ) There are two readily available strategies to explore the error picture of a numerical ODE solver.

   • One proceeds to compute a series of results over a series of doubled step sizes $Lh=(1,2,4,8)·10^{-2}$. Compare these to the expected formula $y_h=y^*+C·h^4+O(h^5)$, that is, apply Richardson extrapolation to compute values for $C$ and $y^*$ and check if they are stable over all the doublings.

   • Another strategy uses MMS (manufactured solutions). Here one would fix some nice function $p$ and consider the equation $$ y''+μ(y^2−1)y'+y=p''+μ(p^2−1)p'+p, ~~~y(0)=p(0),~~y'(0)=p'(0), $$ for instance $p(t)=2\sin(t)$ or with the next perturbation term added. Then $p$ is the exact solution that the numerical solution can be compared against.

   • If you want to actually prove the order, that is, verify the order conditions, read up on Butcher trees, for instance on Butcher's homepage https://www.math.auckland.ac.nz/~butcher/ODE-book-2008/Tutorials/. Or read the original derivation of the low order Runge-Kutta methods in W. Kutta (1901) https://archive.org/details/zeitschriftfrma12runggoog/page/n449