{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0bac6ad1",
   "metadata": {},
   "source": [
    "# Ordinary Differential Equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "dd89c3e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "using SymPy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2e40436",
   "metadata": {},
   "source": [
    "## 1. Separation of Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a10f1cf9",
   "metadata": {},
   "source": [
    "Consider the following example of an initial value problem:\n",
    "\n",
    "$$\n",
    "    \\frac{dy}{dx} = \\frac{x}{y}, \\quad y(3) = 5.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b72f93b",
   "metadata": {},
   "source": [
    "We first formulate this problem with SymPy."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "084800e8",
   "metadata": {},
   "source": [
    "The indepedent variable $x$ has to be declared as a symbol."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7d7c5fae",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$x$"
      ],
      "text/plain": [
       "x"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = Sym(\"x\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec728a75",
   "metadata": {},
   "source": [
    "The dependent variable $y$ is a symbolic function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "40042081",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$y$"
      ],
      "text/plain": [
       "y"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = SymFunction(\"y\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c7d8a84",
   "metadata": {},
   "source": [
    "To define the equation, we must use the ``Eq()``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3b763e3b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\frac{d}{d x} y{\\left(x \\right)} = \\frac{x}{y{\\left(x \\right)}}$"
      ],
      "text/plain": [
       "d           x  \n",
       "--(y(x)) = ----\n",
       "dx         y(x)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ode = Eq(diff(y(x), x), x/y(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "283128e2",
   "metadata": {},
   "source": [
    "We solve the ``ode`` with ``dsolve``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "15ab8505",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\left[\\begin{smallmatrix}y{\\left(x \\right)} = - \\sqrt{C_{1} + x^{2}}\\\\y{\\left(x \\right)} = \\sqrt{C_{1} + x^{2}}\\end{smallmatrix}\\right]$"
      ],
      "text/plain": [
       "2-element Vector{Sym{PyCall.PyObject}}:\n",
       " Eq(y(x), -sqrt(C1 + x^2))\n",
       "  Eq(y(x), sqrt(C1 + x^2))"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sys:1: SymPyDeprecationWarning: \n",
      "\n",
      "non-Expr objects in a Matrix is deprecated. Matrix represents\n",
      "a mathematical matrix. To represent a container of non-numeric\n",
      "entities, Use a list of lists, TableForm, NumPy array, or some\n",
      "other data structure instead.\n",
      "\n",
      "See https://docs.sympy.org/latest/explanation/active-deprecations.html#deprecated-non-expr-in-matrix\n",
      "for details.\n",
      "\n",
      "This has been deprecated since SymPy version 1.9. It\n",
      "will be removed in a future version of SymPy.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "sol = dsolve(ode)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c34fca63",
   "metadata": {},
   "source": [
    "We see there are two solutions, depending on a parameter $C_1$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "644b19cc",
   "metadata": {},
   "source": [
    "The initial condition$y (3) = 5$ has not yet been used.  For that, we need a dictionary to define the initial condition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "cf16cb92",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}\\begin{cases}y{\\left(3 \\right)} & \\text{=>} &5\\\\\\end{cases}\\end{equation*}"
      ],
      "text/plain": [
       "Dict{Sym{PyCall.PyObject}, Int64} with 1 entry:\n",
       "  y(3) => 5"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "initcond = Dict(y(3) => 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "55d7c7e9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$y{\\left(x \\right)} = \\sqrt{x^{2} + 16}$"
      ],
      "text/plain": [
       "          _________\n",
       "         /  2      \n",
       "y(x) = \\/  x  + 16 "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solivp = dsolve(ode, ics=initcond)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fbacf8f",
   "metadata": {},
   "source": [
    "Now we see that we do have a unique solution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b77a566b",
   "metadata": {},
   "source": [
    "## 2. ODEs to for Three Problems"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53325ce2",
   "metadata": {},
   "source": [
    "We apply the scientific method to three problems."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2591d8de",
   "metadata": {},
   "source": [
    "### slowing down"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e549007",
   "metadata": {},
   "source": [
    "*An attack submarine is cruising at 40 knots at a depth of 1000 feet\n",
    "when suddenly the reactor scrams. After 1 minute way has dropped to 30 knots.  How long does the crew have to make repairs before forward motion falls below steerageway of 2 knots?*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1643196",
   "metadata": {},
   "source": [
    "Newton's rule equates the sum of the inertial forces with the sum of the external forces."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d11319a",
   "metadata": {},
   "source": [
    "The inertial force is\n",
    "\n",
    "$$\n",
    "   F = m a,\n",
    "$$\n",
    "\n",
    "where $m$ is the mass and $a$ the acceleration."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02948361",
   "metadata": {},
   "source": [
    "The external force is the water resistance\n",
    "\n",
    "$$\n",
    "    R = -k v^2\n",
    "$$\n",
    "\n",
    "proportional to the square of the velocity, for some constant $k$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "643b8d9f",
   "metadata": {},
   "source": [
    "Note that acceleration is the derivative of the velocity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "bcda53c5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$m \\frac{d}{d t} v{\\left(t \\right)} = - k v^{2}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "  d              2   \n",
       "m*--(v(t)) = -k*v (t)\n",
       "  dt                 "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m, k, t = Sym(\"m, k, t\")\n",
    "v = SymFunction(\"v\")\n",
    "ode = Eq(m*diff(v(t),t), -k*v(t)^2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "467d63af",
   "metadata": {},
   "source": [
    "Let lump the constants $m$ and $k$ into one constant, dividing by $m$ and replacing $k/m$ by $K$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "1f510a78",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\frac{d}{d t} v{\\left(t \\right)}$"
      ],
      "text/plain": [
       "d       \n",
       "--(v(t))\n",
       "dt      "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lefteq = ode.lhs()/m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "75102a3a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$- K v^{2}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "    2   \n",
       "-K*v (t)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "K = Sym(\"K\")\n",
    "righteq = subs(ode.rhs()/m, k/m=>K)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "f3766c64",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\frac{d}{d t} v{\\left(t \\right)} = - K v^{2}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "d              2   \n",
       "--(v(t)) = -K*v (t)\n",
       "dt                 "
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newode = Eq(lefteq, righteq)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80e6f16d",
   "metadata": {},
   "source": [
    "Recall that the *submarine is cruising at 40 knots* so\n",
    "the initial condition is $v(0) = 40$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "6c76de33",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}\\begin{cases}v{\\left(0 \\right)} & \\text{=>} &40\\\\\\end{cases}\\end{equation*}"
      ],
      "text/plain": [
       "Dict{Sym{PyCall.PyObject}, Int64} with 1 entry:\n",
       "  v(0) => 40"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "initcond = Dict(v(0)=>40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "b7c3bcf7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$v{\\left(t \\right)} = \\frac{1}{K t + \\frac{1}{40}}$"
      ],
      "text/plain": [
       "           1     \n",
       "v(t) = ----------\n",
       "       K*t + 1/40"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol = dsolve(newode, ics=initcond)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d532b73a",
   "metadata": {},
   "source": [
    "Recall that *after one minute way has dropped to 30 knots*, so\n",
    "we determine $K$ using $v(1) = 30$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "0688df48",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$30 = \\frac{1}{K + \\frac{1}{40}}$"
      ],
      "text/plain": [
       "        1    \n",
       "30 = --------\n",
       "     K + 1/40"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Kequ = subs(sol, t=>1, v(1) => 30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "d873869c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\left[\\begin{smallmatrix}\\frac{1}{120}\\end{smallmatrix}\\right]$"
      ],
      "text/plain": [
       "1-element Vector{Sym{PyCall.PyObject}}:\n",
       " 1/120"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ksol = solve(Kequ, K)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "12633f09",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$v{\\left(t \\right)} = \\frac{1}{\\frac{t}{120} + \\frac{1}{40}}$"
      ],
      "text/plain": [
       "          1    \n",
       "v(t) = --------\n",
       "        t    1 \n",
       "       --- + --\n",
       "       120   40"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vsol = subs(sol, K=>Ksol[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25950070",
   "metadata": {},
   "source": [
    "The original equation was *how long before the steerageway fall below 2 knots?*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "97359e24",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$2 = \\frac{1}{\\frac{t}{120} + \\frac{1}{40}}$"
      ],
      "text/plain": [
       "       1    \n",
       "2 = --------\n",
       "     t    1 \n",
       "    --- + --\n",
       "    120   40"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vequ = subs(vsol, v(t) => 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "02ad96a6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\left[\\begin{smallmatrix}57\\end{smallmatrix}\\right]$"
      ],
      "text/plain": [
       "1-element Vector{Sym{PyCall.PyObject}}:\n",
       " 57"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve(vequ, t)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e892b38",
   "metadata": {},
   "source": [
    "At 57 minutes the speed has dropped to 2 knots.  So the crew has 56 minutes left."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28dc9492",
   "metadata": {},
   "source": [
    "### cooling off"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0785a49",
   "metadata": {},
   "source": [
    "*A house furnace fails on a cold winter's evening when the outside (ambient) temperature of 20 degrees Fahrenheit. Although initially at 70 degree Fahrenheit, the inside temperature has fallen to 65 degrees Fahrenheit after 1 hour.\n",
    "How long before the inside temperature reaches the damaging temperature of 32 degrees Fahrenheit?*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09d2cca7",
   "metadata": {},
   "source": [
    "Newton's law of cooling states that the rate at which the temperature changes is proportional to the temperature gradient driving out the heat."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fac0415b",
   "metadata": {},
   "source": [
    "Let $T(t)$ be the temperature at time $t$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "b81dfc0d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\frac{d}{d t} T{\\left(t \\right)} = k \\left(T{\\left(t \\right)} - 20\\right)$"
      ],
      "text/plain": [
       "d                       \n",
       "--(T(t)) = k*(T(t) - 20)\n",
       "dt                      "
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T = SymFunction(\"T\")\n",
    "ode = Eq(diff(T(t),t), k*(T(t) - 20))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02e6c466",
   "metadata": {},
   "source": [
    "The right hand side of the equation $T(t) - 20$ represents the difference of the current temperature $T(t)$ with the outside temperature of 20 degrees."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71bea3e9",
   "metadata": {},
   "source": [
    "At $t=0$, the temperature is 70 degrees.  So we formulate the initial condition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "30b9a6ba",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}\\begin{cases}T{\\left(0 \\right)} & \\text{=>} &70\\\\\\end{cases}\\end{equation*}"
      ],
      "text/plain": [
       "Dict{Sym{PyCall.PyObject}, Int64} with 1 entry:\n",
       "  T(0) => 70"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "initcond = Dict(T(0) => 70)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "1df607c9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$T{\\left(t \\right)} = 50 e^{k t} + 20$"
      ],
      "text/plain": [
       "           k*t     \n",
       "T(t) = 50*e    + 20"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol = dsolve(ode, ics=initcond)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "236aa8cb",
   "metadata": {},
   "source": [
    "What is the constant $k$?  Recall that *the inside temperature has fallen to 65 degrees Fahrenheit, after 1 hour.* "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "da6e39ff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$65 = 50 e^{k} + 20$"
      ],
      "text/plain": [
       "         k     \n",
       "65 = 50*e  + 20"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kequ = subs(sol, t=>1, T(1)=>65)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "48df8702",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\left[\\begin{smallmatrix}\\log{\\left(\\frac{9}{10} \\right)}\\end{smallmatrix}\\right]$"
      ],
      "text/plain": [
       "1-element Vector{Sym{PyCall.PyObject}}:\n",
       " log(9/10)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kval = solve(kequ, k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "29d3a9e5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$T{\\left(t \\right)} = 50 e^{t \\log{\\left(\\frac{9}{10} \\right)}} + 20$"
      ],
      "text/plain": [
       "           t*log(9/10)     \n",
       "T(t) = 50*e            + 20"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Tsol = subs(sol, k=>kval[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "893f21c6",
   "metadata": {},
   "source": [
    "What is the original question again?  *How long before the inside temperature reaches the damaging temperature of 32 degrees Fahrenheit?*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "8b010739",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$32 = 50 e^{t \\log{\\left(\\frac{9}{10} \\right)}} + 20$"
      ],
      "text/plain": [
       "         t*log(9/10)     \n",
       "32 = 50*e            + 20"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Tequ = subs(Tsol, T(t)=>32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "f493686f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\left[\\begin{smallmatrix}\\log{\\left(\\left(\\frac{6}{25}\\right)^{\\frac{1}{\\log{\\left(\\frac{9}{10} \\right)}}} \\right)}\\end{smallmatrix}\\right]$"
      ],
      "text/plain": [
       "1-element Vector{Sym{PyCall.PyObject}}:\n",
       " log((6/25)^(1/log(9/10)))"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tval = solve(Tequ, t)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1082dede",
   "metadata": {},
   "source": [
    "To interpret the result, we convert to a 64-bit floating-point approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "c11e5606",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "13.545077553292497"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Float64(tval[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78615fa2",
   "metadata": {},
   "source": [
    "So it takes about 13 and a half hour before it starts freezing inside."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "613e5c39",
   "metadata": {},
   "source": [
    "### population modeling"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b429d45a",
   "metadata": {},
   "source": [
    "*Predict the future population of a developed country.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a71a989",
   "metadata": {},
   "source": [
    "Considering the births, the growth of the population if proportional to its size."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90921abe",
   "metadata": {},
   "source": [
    "Let $P(t)$ be the size of population at time $t$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "f204ddc9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\frac{d}{d t} P{\\left(t \\right)} = k P{\\left(t \\right)}$"
      ],
      "text/plain": [
       "d                \n",
       "--(P(t)) = k*P(t)\n",
       "dt               "
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P = SymFunction(\"P\")\n",
    "ode = Eq(diff(P(t),t), k*P(t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "f12d78cf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$P{\\left(t \\right)} = C_{1} e^{k t}$"
      ],
      "text/plain": [
       "           k*t\n",
       "P(t) = C1*e   "
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dsolve(ode)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31f9ba0e",
   "metadata": {},
   "source": [
    "The exponential growth is not realistic.  Deaths occur."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e9733f4",
   "metadata": {},
   "source": [
    "For a better model, consider the correction term, similar to the second term in a Taylor series.  This second term is proportional to the square of the population size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "4563409d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\frac{d}{d t} P{\\left(t \\right)} = - K P^{2}{\\left(t \\right)} + k P{\\left(t \\right)}$"
      ],
      "text/plain": [
       "d               2            \n",
       "--(P(t)) = - K*P (t) + k*P(t)\n",
       "dt                           "
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "logistics = Eq(diff(P(t),t), k*P(t) - K*(P(t))^2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "615c52a1",
   "metadata": {},
   "source": [
    "This model is called the *logistics equation*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "cdd53693",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$P{\\left(t \\right)} = \\frac{k e^{k \\left(C_{1} + t\\right)}}{K \\left(e^{k \\left(C_{1} + t\\right)} - 1\\right)}$"
      ],
      "text/plain": [
       "             k*(C1 + t)   \n",
       "          k*e             \n",
       "P(t) = -------------------\n",
       "         / k*(C1 + t)    \\\n",
       "       K*\\e           - 1/"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Psol = dsolve(logistics)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "def79ad1",
   "metadata": {},
   "source": [
    "To interpret this model, let us normalize the size of the population at the beginning of time to one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "cb2719b5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$P{\\left(t \\right)} = \\frac{k e^{k \\left(t + \\frac{\\log{\\left(\\frac{K}{K - k} \\right)}}{k}\\right)}}{K \\left(e^{k \\left(t + \\frac{\\log{\\left(\\frac{K}{K - k} \\right)}}{k}\\right)} - 1\\right)}$"
      ],
      "text/plain": [
       "               /       /  K  \\\\   \n",
       "               |    log|-----||   \n",
       "               |       \\K - k/|   \n",
       "             k*|t + ----------|   \n",
       "               \\        k     /   \n",
       "          k*e                     \n",
       "P(t) = ---------------------------\n",
       "         /   /       /  K  \\\\    \\\n",
       "         |   |    log|-----||    |\n",
       "         |   |       \\K - k/|    |\n",
       "         | k*|t + ----------|    |\n",
       "         |   \\        k     /    |\n",
       "       K*\\e                   - 1/"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "initcond = Dict(P(0)=>1)\n",
    "Psol1 = dsolve(logistics, ics=initcond)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93b13055",
   "metadata": {},
   "source": [
    "Fortunately, there is the ``simplify`` we can apply to an expression. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "5de1f908",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\frac{k e^{k t}}{K e^{k t} - K + k}$"
      ],
      "text/plain": [
       "       k*t    \n",
       "    k*e       \n",
       "--------------\n",
       "   k*t        \n",
       "K*e    - K + k"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Psol2 = simplify(Psol1.rhs())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96aef932",
   "metadata": {},
   "source": [
    "Both constants $k$ and $K$ are positive.  Let us divide numerator and denominator by ``exp(k*t)``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "029311ef",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\frac{k}{K + \\left(- K + k\\right) e^{- k t}}$"
      ],
      "text/plain": [
       "        k         \n",
       "------------------\n",
       "              -k*t\n",
       "K + (-K + k)*e    "
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Psol3 = k/(K + (k-K)*exp(-k*t))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51cbe2a0",
   "metadata": {},
   "source": [
    "In the limit, as $t$ goes to $\\infty$, the size of the population reaches $k/K$.  The $k/K$ is called the *carrying capacity* of the society."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.12",
   "language": "julia",
   "name": "julia-1.12"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}