09_Differential_Equations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Differential Equations\n",
    "\n",
    "Let's first revise first-order and second-order differential equations.\n",
    "\n",
    "Second-order differential equations will only be examined for constant coefficients."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First Order Differential Equation:\n",
    "\n",
    "From a differential equation:\n",
    "\n",
    "$$\n",
    "y' + p(x)y = g(x)\n",
    "$$\n",
    "\n",
    "We can solve using the _Integrating Factor_ method.\n",
    "\n",
    "We define $I$ as the integrating factor, equal to:\n",
    "\n",
    "$$\n",
    "I:= e^{\\int p(x) \\ dx}\n",
    "$$\n",
    "\n",
    "We then take $I$ and multiply by it on both sides of our differential equation:\n",
    "\n",
    "$$\n",
    "I y' + I p(x)y = I g(x)\n",
    "$$\n",
    "\n",
    "By taking the derivative of $I$, we can see that:\n",
    "\n",
    "$$\n",
    "I' = p(x) I\n",
    "$$\n",
    "\n",
    "Which if we substitute it into the differential equation, we can see:\n",
    "\n",
    "$$\n",
    "I y' + I'y = I g(x)\n",
    "$$\n",
    "\n",
    "This looks an awful lot like the product rule.\n",
    "\n",
    "We can transform our DE further.\n",
    "\n",
    "$$\n",
    "\\frac{d}{dx} \\left( I y \\right) = I g(x)\n",
    "$$\n",
    "\n",
    "Integrating both sides:\n",
    "\n",
    "$$\n",
    "\\int \\frac{d}{dx} \\left( I y \\right) \\ dx = \\int I g(x) \\ dx\n",
    "$$\n",
    "\n",
    "By the fundamental theorem of calculus, we can simplify the LHS:\n",
    "\n",
    "$$\n",
    "I y = \\int I g(x) \\ dx\n",
    "$$\n",
    "\n",
    "And further:\n",
    "\n",
    "$$\n",
    "y = \\frac{1}{I} \\left( \\int I g(x) \\ dx \\right) \n",
    "\\quad \\equiv \\quad \n",
    "y = \\frac{1}{e^{\\int p(x) \\ dx}} \\left( \\int e^{\\int p(x) \\ dx} g(x) \\ dx \\right)\n",
    "$$\n",
    "\n",
    "Meaning we can (almost) always find the analytic solution to a first order differential equation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Second Order Differential Equations\n",
    "\n",
    "From a differential equation:\n",
    "\n",
    "$$\n",
    "y'' + p(x) y' + q(x) y = g(x)\n",
    "$$\n",
    "\n",
    "We wish to obtain the general solution. \n",
    "\n",
    "The full, general solution is given by:\n",
    "\n",
    "$$\n",
    "y(x) = y_c + y_p\n",
    "$$\n",
    "\n",
    "Where:\n",
    "\n",
    "$$\n",
    "y_c := \\text{ The complementary solution}\n",
    "$$\n",
    "\n",
    "and \n",
    "\n",
    "$$\n",
    "y_p := \\text{ The particular solution}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We generate an auxiliary equation:\n",
    "\n",
    "$$\n",
    "m^2 + p m + q = 0\n",
    "$$\n",
    "\n",
    "A second order polynomial produces (at most) two roots.\n",
    "\n",
    "\n",
    "### Distinct, Real Roots\n",
    "\n",
    "If $m$ is two distinct, real roots, we have:\n",
    "\n",
    "$$\n",
    "y(x) = A e^{m_1 x} + B e^{m_2 x}\n",
    "$$\n",
    "\n",
    "### Repeated, Real Roots\n",
    "\n",
    "If $m$ is two repeated, real roots, we have:\n",
    "\n",
    "$$\n",
    "y(x) = A e^{m x} + B x e^{m x}\n",
    "$$\n",
    "\n",
    "### Complex Roots\n",
    "\n",
    "If $m$ is a pair of complex roots, of the form $p \\pm qi$, we have: \n",
    "\n",
    "$$\n",
    "y(x) = e^{px} \\left( A \\cos (qx) + B \\sin (qx) \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving DEs with `sympy`\n",
    "\n",
    "Setup:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import Symbol, Function, dsolve, diff, init_printing\n",
    "NotImplementedError\n",
    "init_printing()\n",
    "x = Symbol('x')\n",
    "y = Function('y')(x)\n",
    "yp = diff(y)\n",
    "ypp = diff(yp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Say we intend to solve:\n",
    "\n",
    "$$\n",
    "y'' + 5y = 5\n",
    "$$\n",
    "\n",
    "As must be done with the `factor` function we must set the equation equal to zero, e.g.:\n",
    "\n",
    "\n",
    "$$\n",
    "y'' + 5y -5 = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "######################################\n",
    "eq = ypp + 5*y -5 \n",
    "######################################\n",
    "\n",
    "try:\n",
    "    sol = dsolve(eq, y)\n",
    "    display(sol)\n",
    "except:\n",
    "    display(\"Didn't work.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The specific function is `dsolve`, who you can read more about by running the cell below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import dsolve\n",
    "help(dsolve)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 0:\n",
    "\n",
    "Solve the following differential equations using `sympy`:\n",
    "\n",
    "(i) $y' + xy = x^2$\n",
    "\n",
    "(ii) $y'' + 11y' + 30y = 0$\n",
    "\n",
    "(iii) $y'' + 11y' + 30y = 2 \\cos(t) $\n",
    "\n",
    "(iv) $y'' + a^2y = 5 \\quad | \\quad a:= \\text{ const. }$\n",
    "\n",
    "----\n",
    "\n",
    "Below, the equations are not DEs you are likely used to, and are solved with methods you have (probably) not learned yet.\n",
    "\n",
    "For completeness, the methods used to solve these can be _Variation of Parameters_. But, you can solve them here with `sympy`.\n",
    "\n",
    "(v) $y'' - y = \\frac{2}{e^x + 1}$\n",
    "\n",
    "(vi) $y'' + y = \\tan(x) \\sec(x)$\n",
    "\n",
    "## Answer 0:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `try` and `except` blocks.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the following code:\n",
    "\n",
    "### Example 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import Symbol, Function, dsolve, diff, init_printing, tan, sec, csc, E as e\n",
    "init_printing()\n",
    "x = Symbol('x')\n",
    "y = Function('y')(x)\n",
    "yp = diff(y)\n",
    "ypp = diff(yp)\n",
    "\n",
    "######################################\n",
    "eq = (x**3)*ypp + x*yp -y \n",
    "######################################\n",
    "\n",
    "dsolve(eq, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see the code fails horribly with something called a `NotImplementedError`. If solving this DE were one of a sequence of steps, it would cause our code to exit and the program to halt -- not particularly desirable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 1:\n",
    "\n",
    "Looking back at the code before question 0, we can see a use of `try` and `except`. \n",
    "\n",
    "(i) What do you think the use of these keywords does? \n",
    "\n",
    "_(If this perplexes you, try commenting out the lines and changing the indents to see what happens without it)._\n",
    "\n",
    "(ii) Attempt to re-implement the code in `Example 1` below with a `try`-`except` block.\n",
    "\n",
    "## Answer 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Systems of Differential Equations\n",
    "\n",
    "Consider the system:\n",
    "\n",
    "$$\n",
    "\\frac{dx_1}{dt} = -2 x_1 + 5 x_2 \n",
    "\\quad;\\quad\n",
    "\\frac{dx_2}{dt} = x_1 + 2 x_2 \n",
    "$$\n",
    "\n",
    "We can form this as a matrix:\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "    x_1' \\\\\n",
    "    x_2' \\\\\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "    -2 & 5 \\\\\n",
    "    1 & 2 \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "    x_1 \\\\\n",
    "    x_2\\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "And obtain the solution as:\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "    x_1 \\\\\n",
    "    x_2\\\\\n",
    "\\end{bmatrix}\n",
    "=\n",
    "C_1 e^{\\lambda_1 t} \\mathbf{v}_1 + C_2 e^{\\lambda_2 t} \\mathbf{v}_2\n",
    "$$\n",
    "\n",
    "Where $\\lambda_n$ is the $n$-th eigenvalue and $\\mathbf{v}_n$ is the $n$-th eigen-vector."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's generate the matrix to get started:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import Matrix, Symbol, init_printing\n",
    "init_printing()\n",
    "\n",
    "x1 = Symbol('x_1')\n",
    "x2 = Symbol('x_2')\n",
    "\n",
    "A = Matrix([\n",
    "    [-2,5],\n",
    "    [1,2]\n",
    "])\n",
    "display(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution to the system is heavily dependent on whether we can diagonalise the matrix.\n",
    "\n",
    "i.e.\n",
    "\n",
    "$$\n",
    "\\mathbf{A} = \\mathbf{PDP^{-1}}\n",
    "$$\n",
    "\n",
    "Where $\\mathbf{D}$ is a diagonal matrix.\n",
    "\n",
    "We can easily diagonalise with a simple `method` called `diagonalize()` (note the american-english spelling)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P,D = A.diagonalize()\n",
    "\n",
    "display('P')\n",
    "display(P)\n",
    "display('D')\n",
    "display(D)\n",
    "display('P^-1')\n",
    "display(P**-1)\n",
    "\n",
    "# Are they the same?\n",
    "display('A == P*D*P^-1  ?')\n",
    "if A == P*D*P**-1:\n",
    "    display(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that $\\mathbf{D}$ is in fact a diagonal matrix of the eigenvalues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.eigenvals()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What of the eigenvectors?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.eigenvects()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It appears that $\\mathbf{P}$ is simply our eigenvectors collated into one matrix.\n",
    "\n",
    "So by the assertion at the beginning of this section, the solution must be:\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "    x_1 \\\\\n",
    "    x_2\\\\\n",
    "\\end{bmatrix}\n",
    "=\n",
    "C_1 e^{\\lambda_1 t} \\mathbf{v}_1 + C_2 e^{\\lambda_2 t} \\mathbf{v}_2\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "    x_1 \\\\\n",
    "    x_2\\\\\n",
    "\\end{bmatrix}\n",
    "=\n",
    "C_1 e^{-3 t} \\begin{bmatrix} -5 \\\\ 1 \\end{bmatrix} + C_2 e^{3 t} \\begin{bmatrix} 1 \\\\ 1 \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 2\n",
    "\n",
    "What are the solutions to the following systems?\n",
    "\n",
    "(i) \n",
    "\n",
    "$$\n",
    "\\frac{dx}{dt} = x - 8y \n",
    "\\quad ; \\quad\n",
    "\\frac{dy}{dt} = 2x + y\n",
    "$$\n",
    "\n",
    "(ii) \n",
    "\n",
    "$$\n",
    "\\frac{dx}{dt} = x + 2y - 3z\n",
    "\\quad ; \\quad\n",
    "\\frac{dy}{dt} = -5x + y - 4z\n",
    "\\quad ; \\quad\n",
    "\\frac{dz}{dt} = -2y + 4z\n",
    "$$\n",
    "\n",
    "## Answer 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}