06_ODE_Ex.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. Ordinary Differential Equations: Exercise\n",
    "\n",
    "Name: \n",
    "\n",
    "Date: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from scipy.integrate import odeint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Linear ODEs\n",
    "Like the exponential of a real number $x$ is given by\n",
    "\n",
    "$$ e^{x} = 1 + x + \\frac{1}{2} x^2 + \\frac{1}{6} x^3 + ...\n",
    "  = \\sum_{k=0}^{\\infty} \\frac{1}{k!} x^k, $$\n",
    "\n",
    "the *exponential of a matrix* $X$ is defined as\n",
    "\n",
    "$$ e^{X} = I + X + \\frac{1}{2} X^2 + \\frac{1}{6} X^3 + ...\n",
    "  = \\sum_{k=0}^{\\infty} \\frac{1}{k!} X^k. $$\n",
    "\n",
    "For one dimensional linear ODE\n",
    "\n",
    "$$ \\frac{dy}{dt} = a y $$\n",
    "\n",
    "the solution is given by\n",
    "\n",
    "$$ y(t) = e^{at} y(0), $$\n",
    "\n",
    "where $y(0)$ is the initial state. \n",
    "\n",
    "For an $n$ dimensional linear ODE\n",
    "\n",
    "$$ \\frac{dy}{dt} = A y $$\n",
    "\n",
    "where $A$ is an $n\\times n$ matrix, the solution is given by the matrix exponential\n",
    "\n",
    "$$ y(t) = e^{At} y(0), $$\n",
    "\n",
    "where $y(0)$ is an $n$-dimensional initial state.\n",
    "\n",
    "* Verify this by expanding $e^{At}$ accordint to the definition and differentiating each term by $t$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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.\n",
    "\n",
    "Let us visualize solutions for different eigenvalues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def linear(y, t, A):\n",
    "    \"\"\"Linear dynamcal system dy/dt = Ay\n",
    "    y: n-dimensional state vector\n",
    "    t: time (not used, for compatibility with odeint())\n",
    "    A: n*n matrix\"\"\"\n",
    "    # y is an array (row vector), A is a matrix\n",
    "    return A@y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def linear2D(A, yinit=np.array([[1,0],[0,1],[-1,0],[0,-1]]), t=np.arange(0, 5, 0.1)):\n",
    "    \"\"\"Visualizing linear 2D dynamical system\"\"\"\n",
    "    for y0 in yinit:\n",
    "        y = odeint(linear, y0, t, args=(A,))\n",
    "        plt.plot(y[0,0], y[0,1], 'o')   # starting point\n",
    "        plt.plot(y[:,0], y[:,1], '+-')  # trajectory\n",
    "    plt.axis('equal')\n",
    "    return np.linalg.eig(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "1) Real eigenvalues $\\lambda_1 > \\lambda_2 > 0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[2, 0], [0, 1]])  # modify this!\n",
    "linear2D(A, t=np.arange(0, 1, 0.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "2) Real eigenvalues $\\lambda_1 > 0 > \\lambda_2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "3) Real eigenvalues $0 > \\lambda_1 > \\lambda_2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4) Complex eigenvalues $\\lambda_1=a+ib$ and $\\lambda_2=a-ib$ with $a>0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5) Complex eigenvalues $\\lambda_1=a+ib$ and $\\lambda_2=a-ib$ with $a<0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "c.f. For a 2 by 2 matrix\n",
    "\n",
    "$$ A = \\pmatrix{a & b \\\\ c & d}, $$\n",
    "\n",
    "we can analytically derive the eivenvalues from \n",
    "\n",
    "$$ \\det (A - \\lambda I) = (a-\\lambda)(d-\\lambda) - bc = 0 $$\n",
    "\n",
    "as\n",
    "\n",
    "$$ \\lambda = \\frac{a+d}{2} \\pm \\sqrt{\\frac{(a-d)^2}{4}+ bc}. $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Nonlinear ODEs\n",
    "1) Implement a nonlinear system, such as a pendulum with friction $\\mu$:\n",
    "\n",
    "$$ \\frac{d\\theta}{dt} = \\omega $$\n",
    "\n",
    "$$ ml^2 \\frac{d\\omega}{dt} = - \\mu \\omega - mgl \\sin \\theta $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pendulum(y, t, m=1, l=1, mu=1, g=9.8):\n",
    "    \"\"\"pendulum dynamics\n",
    "    m:mass, l:length, mu:damping, g:gravity\"\"\"\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2) Run a simulation by `odeint()` and show the trajectory as (t, y(t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3) Show the trajectory in the 2D state space $(\\theta, \\omega)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Option) Implement a nonlinear system with time-dependent input, such as a forced pendulum:\n",
    "\n",
    "$$ \\frac{d\\theta}{dt} = \\omega $$\n",
    "\n",
    "$$ ml^2 \\frac{d\\omega}{dt} = - \\mu \\omega - mgl \\sin\\theta + a\\sin bt$$\n",
    "\n",
    "and see how the behavior changes with the input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Bifurcation\n",
    "FitzHugh-Nagumo model is an extension of Van der Pol model to approximate spiking behaviors of neurons.\n",
    "\n",
    "$$ \\frac{dv}{dt} = v - \\frac{v^3}{3} - w + I $$\n",
    "\n",
    "$$ \\frac{dw}{dt} = \\phi (v + a - bw) $$\n",
    "\n",
    "1) Implement a function and see how the behaviors at different input current $I$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fhn(y, t, I=0, a=0.7, b=0.8, phi=0.08):\n",
    "    \"\"\"FitzHugh-Nagumo model\"\"\"\n",
    "    v, w = y\n",
    "\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y0 = np.array([0, 0])\n",
    "t = np.arange(0, 100, 0.1)\n",
    "y = odeint(fhn, ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(y[:,0], y[:,1], '+-')  # phase plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2) Draw a *bifurcation diagram* showing the max-min of $v$ for different values of $I$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}