Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Van der pol oscillator
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment