Home Statistics and Visualization The Endeavour Solving Van der Pol equation with ivp_solve

Solving Van der Pol equation with ivp_solve

The Endeavour 2019-12-22

Van der Pol’s differential equation is

{d^2x \over dt^2}-\mu(1-x^2){dx \over dt}+x= 0

The equation describes a system with nonlinear damping, the degree of damping given by μ. If μ = 0 the system is linear and undamped, but for positive μ the system is nonlinear and damped. We will plot the phase portrait for the solution to Van der Pol’s equation in Python using SciPy’s new ODE solver ivp_solve.

The function ivp_solve does not solve second-order systems of equations directly. It solves systems of first-order equations, but a second-order differential equation can be recast as a pair of first-order equations by introducing the first derivative as a new variable.

\begin{align*} {dx \over dt} &= y \\ {dy \over dt}&= \mu(1-x^2)y -x \\ \end{align*}

Since y is the derivative of x, the phase portrait is just the plot of (x, y).

Phase portait of Van der Pol oscillator

If μ = 0, we have a simple harmonic oscillator and the phase portrait is simply a circle. For larger values of μ the solutions enter limiting cycles, but the cycles are more complicated than just circles.

Here’s the Python code that made the plot.

from scipy import linspacefrom scipy.integrate import solve_ivpimport matplotlib.pyplot as pltdef vdp(t, z):    x, y = z    return [y, mu*(1 - x**2)*y - x]a, b = 0, 10mus = [0, 1, 2]styles = ["-", "--", ":"]t = linspace(a, b, 1000)for mu, style in zip(mus, styles):    sol = solve_ivp(vdp, [a, b], [1, 0], t_eval=t)    plt.plot(sol.y[0], sol.y[1], style)  # make a little extra horizontal room for legendplt.xlim([-3,3])    plt.legend([f"$\mu={m}$" for m in mus])plt.axes().set_aspect(1)

Related posts