The Duffing oscillator

The Duffing oscillator is the second-order differential equation $\ddot{x}=2x-x^3$. The function solve_ivp() does not deal directly with second-order equations, but there is a standard way to reformulate them as first-order equations with extra variables. For the Duffing oscillator, we simply define $y$ to be $\dot{x}$ so that $\dot{y}=\ddot{x}=2x-x^3$ and $\dot{x}=y$. We now have a first order equation for the vector $u=\begin{bmatrix} x \\ y \end{bmatrix}$, namely $$ \dot{u} = \begin{bmatrix} \dot{x} \\ \dot{y} \end{bmatrix} = \begin{bmatrix} y \\ 2x-x^3 \end{bmatrix} = \begin{bmatrix} u_1 \\ 2u_0-u_0^3 \end{bmatrix} $$

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

We plot the solutions starting at various different initial points, some on the positive $x$-axis, some on the negative $x$-axis, and some on the positive $y$-axis. There are three possible constant solutions: one where $x=\dot{x}=0$ for all $t$, one where $x=\sqrt{2}$ and $\dot{x}=0$ for all $t$, and one where $x=-\sqrt{2}$ and $\dot{x}=0$ for all $t$. These are marked with black dots.

def du_dt(t, u):
    return np.array([u[1], 2*u[0] - u[0]**3])

fig, ax = plt.subplots(figsize = (8, 8))
ax.set_aspect('equal')

t = np.linspace(0, 24, 4000)
ics = ([[ x0, 0] for x0 in np.arange(0.1, 1.5,0.1)] + 
       [[-x0, 0] for x0 in np.arange(0.1, 1.5,0.1)] + 
       [[ 0, y0] for y0 in np.arange(0.2, 1.5,0.1)])
for xy0 in ics:
    sol = solve_ivp(du_dt, [t[0], t[-1]], xy0, t_eval = t, method = 'DOP853', atol = 1e-9)
    ax.plot(sol.y[0], sol.y[1])
ax.plot([-np.sqrt(2), 0, np.sqrt(2)], [0, 0, 0], 'ko')

[<matplotlib.lines.Line2D at 0x23439f7b350>]