ODEs¶

In [8]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
#%matplotlib osx
In [9]:
X, Y = np.meshgrid( np.linspace(-1, 1., 21), np.linspace(-1,1., 21))

theta = 0.9 * 2 * np.pi
alpha = 0.4
U = -theta * Y - alpha*X
V = theta * X - alpha*Y

plt.figure(1)
plt.clf()
Q = plt.quiver(X,Y,U,V,units='xy', width=0.005, headwidth=4)
plt.axis('equal');

The function below implements the decaying harmonic oscillator system of ODEs \begin{align} x' &= -\theta y - \alpha x \\ y' &= \theta x - \alpha y \end{align}

In [10]:
# This function evaluates the derivatives of all the dynamic variables.
# ie. it evaluates the RHS of the system of ODEs.

def fcn(t, x):
    return [-theta*x[1] - alpha*x[0] , theta*x[0] - alpha*x[1]]
In [11]:
plt.figure(1)
plt.clf()
Q = plt.quiver(X,Y,U,V,units='xy', width=0.005, headwidth=4)
plt.axis('equal');

tspan = [0, 3.]
y0 = [0, 0.9]
sol = solve_ivp(fcn, tspan, y0, max_step=0.01)

plt.plot(sol.y[0], sol.y[1], 'r', lw=3);
plt.plot(sol.y[0,0], sol.y[1,0], 'ro');
In [ ]: