6. Ordinary Differential Equations: Exercise

6. Ordinary Differential Equations: Exercise#

Name:

Date:

(Please submit this .ipynb file with your name and its PDF copy.)

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, 0], [0, 1]])  # modify this!
linear2D(A, t=np.arange(0, 1, 0.1))

EigResult(eigenvalues=array([2., 1.]), eigenvectors=array([[1., 0.],
       [0., 1.]]))
_images/9ecb93c1cde06684aef4b8b569e70bca105b82486ff6fd916453b129c9cde1c3.png

  1. Real eigenvalues \(\lambda_1 > 0 > \lambda_2\)

  1. Real eigenvalues \(0 > \lambda_1 > \lambda_2\)

  1. Complex eigenvalues \(\lambda_1=a+ib\) and \(\lambda_2=a-ib\) with \(a>0\)

  1. Complex eigenvalues \(\lambda_1=a+ib\) and \(\lambda_2=a-ib\) with \(a<0\)

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=1, g=9.8):
    """pendulum dynamics
    m:mass, l:length, mu:damping, g:gravity"""

  1. Run a simulation by odeint() and show the trajectory as (t, y(t))

  1. Show the trajectory in the 2D state space \((\theta, \omega)\)

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

    
    

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

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[7], line 3
      1 y0 = np.array([0, 0])
      2 t = np.arange(0, 100, 0.1)
----> 3 y = odeint(fhn, y0, t, args=(0.5,))
      4 plt.plot(t, y, '.-')

File /opt/anaconda3/lib/python3.12/site-packages/scipy/integrate/_odepack_py.py:243, in odeint(func, y0, t, args, Dfun, col_deriv, full_output, ml, mu, rtol, atol, tcrit, h0, hmax, hmin, ixpr, mxstep, mxhnil, mxordn, mxords, printmessg, tfirst)
    241 t = copy(t)
    242 y0 = copy(y0)
--> 243 output = _odepack.odeint(func, y0, t, args, Dfun, col_deriv, ml, mu,
    244                          full_output, rtol, atol, tcrit, h0, hmax, hmin,
    245                          ixpr, mxstep, mxhnil, mxordn, mxords,
    246                          int(bool(tfirst)))
    247 if output[-1] < 0:
    248     warning_msg = (f"{_msgs[output[-1]]} Run with full_output = 1 to "
    249                    f"get quantitative information.")

RuntimeError: The size of the array returned by func (1) does not match the size of y0 (2).

plt.plot(y[:,0], y[:,1], '+-')  # phase plot

  1. Draw a bifurcation diagram showing the max-min of \(v\) for different values of \(I\).

Contents