6. Ordinary Differential Equations: Exercise Solutions

6. Ordinary Differential Equations: Exercise Solutions#

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.integrate import odeint

1. Linear ODEs#

Like the exponential of a real number \(x\) is given by $\( e^{x} = 1 + x + \frac{1}{2} x^2 + \frac{1}{6} x^3 + ... = \sum_{k=0}^{\infty} \frac{1}{k!} x^k, \)\( the *exponential of a matrix* \)X\( is defined as \)\( e^{X} = I + X + \frac{1}{2} X^2 + \frac{1}{6} X^3 + ... = \sum_{k=0}^{\infty} \frac{1}{k!} X^k. \)$

For one dimensional linear ODE $\( \frac{dy}{dt} = a y \)\( the solution is given by \)\( y(t) = e^{at} y(0), \)\( where \)y(0)$ is the initial state.

For an \(n\) dimensional linear ODE $\( \frac{dy}{dt} = A y \)\( where \)A\( is an \)n\times n\( matrix, the solution is given by the matrix exponential \)\( y(t) = e^{At} y(0), \)\( where \)y(0)\( is an \)n$-dimensional initial state.

  • Verify this by expanding \(e^{At}\) accordint to the definition and differentiating each term by \(t\).

The behavior of the matrix exponentioal \(e^{At}\) depends on the eivenvalues of \(A\); whether the eigenvalues are real or complex, and whether the real part is positive or negative.

Let us visualize solutions for different eigenvalues.

def linear(y, t, A):
    """Linear dynamcal system dy/dt = Ay
    y: n-dimensional state vector
    t: time (not used, for compatibility with odeint())
    A: n*n matrix"""
    # y is an array (row vector), A is a matrix
    return A@y
def linear2D(A, yinit=np.array([[1,0],[0,1],[-1,0],[0,-1]]), t=np.arange(0, 5, 0.1)):
    """Visualizing linear 2D dynamical system"""
    for y0 in yinit:
        y = odeint(linear, y0, t, args=(A,))
        plt.plot(y[0,0], y[0,1], 'o')   # starting point
        plt.plot(y[:,0], y[:,1], '+-')  # trajectory
    plt.axis('equal')
    return np.linalg.eig(A)
  1. Real eigenvalues \(\lambda_1 > \lambda_2 > 0\)

A = np.array([[2, 1], [0.5, 1]])  # modift this!
linear2D(A, t=np.arange(0, 1, 0.1))
EigResult(eigenvalues=array([2.3660254, 0.6339746]), eigenvectors=array([[ 0.9390708 , -0.59069049],
       [ 0.34372377,  0.80689822]]))
_images/bef76d5357f3fe3ff732c3ba512b763a05a9bc56c8c9fdecc0ff7b5393c95c77.png
  1. Real eigenvalues \(\lambda_1 > 0 > \lambda_2\)

A = np.array([[2, 1], [0.5, -1]])
linear2D(A, t=np.arange(0, 1, 0.1))
EigResult(eigenvalues=array([ 2.1583124, -1.1583124]), eigenvectors=array([[ 0.98769933, -0.30185542],
       [ 0.15636505,  0.95335371]]))
_images/2c79b9d1e6881bff2848c3439f16f43f7837eff91ea8fc2b79e4f2f971f768f1.png
  1. Real eigenvalues \(0 > \lambda_1 > \lambda_2\)

A = np.array([[-2, 1], [0.5, -1]])
linear2D(A, t=np.arange(0, 1, 0.1))
EigResult(eigenvalues=array([-2.3660254, -0.6339746]), eigenvectors=array([[-0.9390708 , -0.59069049],
       [ 0.34372377, -0.80689822]]))
_images/5c82142b120aea27b96f533b67010c60dd7f28a31d70eee4297d42f9b8b4f1c9.png
  1. Complex eigenvalues \(\lambda_1=a+ib\) and \(\lambda_2=a-ib\) with \(a>0\)

A = np.array([[2, 2], [-1, 1]])
linear2D(A, t=np.arange(0, 1, 0.1))
EigResult(eigenvalues=array([1.5+1.32287566j, 1.5-1.32287566j]), eigenvectors=array([[ 0.81649658+0.j        ,  0.81649658-0.j        ],
       [-0.20412415+0.54006172j, -0.20412415-0.54006172j]]))
_images/73b519fd95584b28e62f42a83d642e1f9f53180f226678897691a6bab15ab6d2.png
  1. Complex eigenvalues \(\lambda_1=a+ib\) and \(\lambda_2=a-ib\) with \(a<0\)

A = np.array([[-2, 1], [-1, -1]])
linear2D(A, t=np.arange(0, 1, 0.1))
EigResult(eigenvalues=array([-1.5+0.8660254j, -1.5-0.8660254j]), eigenvectors=array([[0.70710678+0.j        , 0.70710678-0.j        ],
       [0.35355339+0.61237244j, 0.35355339-0.61237244j]]))
_images/3f65ad3f9ed03d4ca1e0a218b159f44ff18d250ce703a95f3de7872d70dd699a.png

c.f. For a 2 by 2 matrix $\( A = \pmatrix{a & b \\ c & d}, \)\( we can analytically derive the eivenvalues from \)\( \det (A - \lambda I) = (a-\lambda)(d-\lambda) - bc = 0 \)\( as \)\( \lambda = \frac{a+d}{2} \pm \sqrt{\frac{(a-d)^2}{4}+ bc}. \)$

2. Nonlinear ODEs#

  1. Implement a nonlinear system, such as a pendulum with friction \(\mu\): $\( \frac{d\theta}{dt} = \omega \)\( \)\( ml^2 \frac{d\omega}{dt} = - \mu \omega - mgl \sin \theta \)$

def pendulum(y, t, m=1, l=1, mu=0.1, g=9.8):
    th, om = y
    dthdt = om
    domdt = -(mu*om + m*g*l*np.sin(th))/(m*l*l)
    return [dthdt, domdt]
  1. Run a simulation by odeint() and show the trajectory as (t, y(t))

y0 = np.array([3, 0])
t = np.arange(0, 20, 0.1)
y = odeint(pendulum, y0, t, args=(1, 1, 0.2, 9.8))
plt.plot(t, y)  # trajectory
plt.xlabel('t'); plt.ylabel('$\\theta$, $\omega$');
<>:5: SyntaxWarning: invalid escape sequence '\o'
<>:5: SyntaxWarning: invalid escape sequence '\o'
/var/folders/dd/3kwwsm055d1cc1w134st1f580000gn/T/ipykernel_62061/284280705.py:5: SyntaxWarning: invalid escape sequence '\o'
  plt.xlabel('t'); plt.ylabel('$\\theta$, $\omega$');
_images/c160a12c392fe8a990db592ae581212f2a647e03d092b25d87fe725fee223158.png
  1. Show the trajectory in the 2D state space \((\theta, \omega)\)

plt.plot(y[:,0], y[:,1]);  # phase plot
plt.xlabel('$\\theta$'); plt.ylabel('$\omega$');
<>:2: SyntaxWarning: invalid escape sequence '\o'
<>:2: SyntaxWarning: invalid escape sequence '\o'
/var/folders/dd/3kwwsm055d1cc1w134st1f580000gn/T/ipykernel_62061/2133674605.py:2: SyntaxWarning: invalid escape sequence '\o'
  plt.xlabel('$\\theta$'); plt.ylabel('$\omega$');
_images/de4446e1baaa908dcfd2054690b7a5cfa51e7fc58ceabf9c3dcdd8ea2be8e92e.png

Option) Implement a nonlinear system with time-dependent input, such as a forced pendulum: $\( \frac{d\theta}{dt} = \omega \)\( \)\( ml^2 \frac{d\omega}{dt} = - \mu \omega - mgl \sin\theta + a\sin bt\)$ and see how the behavior changes with the input.

3. Bifurcation#

FitzHugh-Nagumo model is an extension of Van der Pol model to approximate spiking behaviors of neurons. $\( \frac{dv}{dt} = v - \frac{v^3}{3} - w + I \)\( \)\( \frac{dw}{dt} = \phi (v + a - bw) \)$

  1. Implement a function and see how the behaviors at different input current \(I\).

def fhn(y, t, I=0, a=0.7, b=0.8, phi=0.08):
    """FitzHugh-Nagumo model"""
    v, w = y
    dvdt = v - v**3/3 - w + I
    dwdt = phi*(v + a -b*w)
    return np.array([dvdt, dwdt])
y0 = np.array([0, 0])
t = np.arange(0, 100, 0.1)
y = odeint(fhn, y0, t, args=(0.5,))
plt.plot(t, y, '.-')  # trajectory
plt.xlabel('t'); plt.ylabel('v, w');
_images/14c739fc7958ab38eba7af17a92ec8bd20f7ea8382081d96f54a1d5298a96dc6.png
plt.plot(y[:,0], y[:,1], '+-')  # phase plot
plt.xlabel('v'); plt.ylabel('w');
_images/307839e68077c01944d9b4993c521757d67a97d65ccc1e4f0545472cbeeb9575.png
  1. Draw a bifurcation diagram showing the max-min of \(v\) for different values of \(I\).

# bifurcation diagram
t = np.arange(0, 200, 0.1)
for I in np.arange(0., 2., 0.05):
    y = odeint(fhn, y0, t, args=(I,))  # change I
    plt.plot(I, min(y[1000:,0]), "b.")
    plt.plot(I, max(y[1000:,0]), "r.")
plt.xlabel('I'); plt.ylabel('v');
_images/af6210ea7feb487e7dc621e1bb4911b875bfee47edf994f27169701295c710f1.png