Van der pol oscillator

van_der_pol.py

import scipy as sc
import numpy as np
from scipy import integrate
from matplotlib import pyplot as plt
from matplotlib.ticker import (MultipleLocator, FormatStrFormatter, AutoMinorLocator, AutoLocator)


def osz(y, t, gamma1, gamma2, omega0):
    x, p = y
    dydt = [p, -(-gamma1 + gamma2 * x**2)*p - omega0**2 * x]
    return dydt

def main():
    t = np.arange(0, 100, 0.01)
    sol1 = sc.integrate.odeint(osz, y0=[0, 1], t=t, args=(1, 1, 1))
    sol2 = sc.integrate.odeint(osz, y0=[0, 1], t=t, args=(1, 1, 5))
    sol3 = sc.integrate.odeint(osz, y0=[0, 1], t=t, args=(1, 1, 1/5))

    fig, ax = plt.subplots(1, 1, figsize=(10, 7), dpi=100)

    ax.xaxis.set_minor_locator(AutoMinorLocator())
    ax.xaxis.set_major_locator(AutoLocator())
    ax.yaxis.set_minor_locator(AutoMinorLocator())
    ax.yaxis.set_major_locator(AutoLocator())

    ax.tick_params(which='both', width=0.8)
    ax.tick_params(which='major', length=7)
    ax.tick_params(which='minor', length=4, grid_linestyle="--", grid_color="lightgray")

    ax.grid(True, which="both")

    ax.plot(sol1[:, 0], sol1[:, 1], label="$\gamma_1 / \omega_0 = 1$", linewidth=1)
    ax.plot(sol2[:, 0], sol2[:, 1], label="$\gamma_1 / \omega_0 = 1/5$", linewidth=1)
    ax.plot(sol3[:, 0], sol3[:, 1], label="$\gamma_1 / \omega_0 = 5$", linewidth=1)

    ax.set_title("Van der Pol, $\gamma_1 = \gamma_2 = 1$")

    ax.legend()

    ax.set_xlabel("$x(t)$")
    ax.set_ylabel("$p(t)$")

    plt.show()



if __name__ == "__main__":
    main()