chapters/chap09.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a8b892c8",
   "metadata": {},
   "source": [
    "Printed and electronic copies of *Modeling and Simulation in Python* are available from [No Starch Press](https://nostarch.com/modeling-and-simulation-python) and [Bookshop.org](https://bookshop.org/p/books/modeling-and-simulation-in-python-allen-b-downey/17836697?ean=9781718502161) and [Amazon](https://amzn.to/3y9UxNb)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "asian-customer",
   "metadata": {},
   "source": [
    "# Analysis of Population Growth"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "imported-table",
   "metadata": {
    "tags": []
   },
   "source": [
    "*Modeling and Simulation in Python*\n",
    "\n",
    "Copyright 2021 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "electoral-turkey",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# install Pint if necessary\n",
    "\n",
    "try:\n",
    "    import pint\n",
    "except ImportError:\n",
    "    !pip install pint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "formal-context",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# download modsim.py if necessary\n",
    "\n",
    "from os.path import basename, exists\n",
    "\n",
    "def download(url):\n",
    "    filename = basename(url)\n",
    "    if not exists(filename):\n",
    "        from urllib.request import urlretrieve\n",
    "        local, _ = urlretrieve(url, filename)\n",
    "        print('Downloaded ' + local)\n",
    "    \n",
    "download('https://raw.githubusercontent.com/AllenDowney/' +\n",
    "         'ModSimPy/master/modsim.py')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "progressive-typing",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import functions from modsim\n",
    "\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "biblical-murder",
   "metadata": {},
   "source": [
    "In this chapter we'll express the models from previous chapters as\n",
    "difference equations and differential equations, solve the equations,\n",
    "and derive the functional forms of the solutions. And I'll present some thoughts about the complementary roles of mathematical analysis and simulation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "hybrid-retrieval",
   "metadata": {},
   "source": [
    "## Difference Equations\n",
    "\n",
    "The population models in the previous chapter and this one are simple\n",
    "enough that we didn't really need to run simulations. We could have\n",
    "solved them mathematically. For example, we wrote the constant growth\n",
    "model like this:\n",
    "\n",
    "```\n",
    "results[t+1] = results[t] + annual_growth\n",
    "```\n",
    "\n",
    "In mathematical notation, we would write the same model like this:\n",
    "\n",
    "$$x_{n+1} = x_n + c$$ \n",
    "\n",
    "where $x_n$ is the population during year $n$, \n",
    "$x_{n+1}$ is the population during year $n+1$,\n",
    "and $c$ is constant annual growth.\n",
    "This way of representing the model, where future population is a function of current population, is a *difference equation*; see\n",
    "<https://en.wikipedia.org/wiki/Linear_difference_equation>.\n",
    "\n",
    "For a given value of $n$, sometimes it is possible to compute $x_n$ directly; that\n",
    "is, without computing the intervening values from $x_1$ through $x_{n-1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "specified-mexican",
   "metadata": {},
   "source": [
    "In the case of constant growth we can see that $x_1 = x_0 + c$, and\n",
    "$x_2 = x_1 + c$. Combining these, we get $x_2 = x_0 + 2c$, then\n",
    "$x_3 = x_0 + 3c$, and we can see that in general\n",
    "\n",
    "$$x_n = x_0 + nc$$ \n",
    "\n",
    "So if we want to know $x_{100}$ and we don't care about the other values, we can compute it with one multiplication and one addition."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sufficient-tampa",
   "metadata": {},
   "source": [
    "We can also write the proportional model as a difference equation:\n",
    "\n",
    "$$x_{n+1} = x_n + \\alpha x_n$$ \n",
    "\n",
    "Or more conventionally as:\n",
    "\n",
    "$$x_{n+1} = x_n (1 + \\alpha)$$ \n",
    "\n",
    "Now we can see that $x_1 = x_0 (1 + \\alpha)$, and $x_2 = x_0 (1 + \\alpha)^2$, and in general\n",
    "\n",
    "$$x_n = x_0 (1 + \\alpha)^n$$\n",
    "\n",
    "A sequence with this functional form is called a *geometric progression*;\n",
    "see <http://modsimpy.com/geom>. When $\\alpha$ is positive, the factor\n",
    "$1+\\alpha$ is greater than 1, so the elements of the sequence grow\n",
    "without bound."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "motivated-consumer",
   "metadata": {},
   "source": [
    "Finally, we can write the quadratic model like this:\n",
    "\n",
    "$$x_{n+1} = x_n + \\alpha x_n + \\beta x_n^2$$ \n",
    "\n",
    "or with the more conventional parameterization like this:\n",
    "\n",
    "$$x_{n+1} = x_n + r x_n (1 - x_n / K)$$ \n",
    "\n",
    "There is no analytic solution to this equation, but we can approximate it with a differential equation and solve that, which is what we'll do in the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "decreased-location",
   "metadata": {},
   "source": [
    "## Differential Equations\n",
    "\n",
    "Starting again with the constant growth model \n",
    "\n",
    "$$x_{n+1} = x_n + c$$ \n",
    "\n",
    "If we define $\\Delta x$ to be the change in $x$ from one time step to the next, we can write: \n",
    "\n",
    "$$\\Delta x = x_{n+1} - x_n = c$$ \n",
    "\n",
    "If we define\n",
    "$\\Delta t$ to be the time step, which is one year in the example, we can\n",
    "write the rate of change per unit of time like this:\n",
    "\n",
    "$$\\frac{\\Delta x}{\\Delta t} = c$$ \n",
    "\n",
    "This is a *discrete* model, which\n",
    "means time is only defined at integer values of $n$ and not in between.\n",
    "But in reality, people are born and die all the time, not once a year,\n",
    "so it might be more realistic to use a *continuous* model, which means\n",
    "time is defined at all values of $t$, not just integers. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "latest-impression",
   "metadata": {},
   "source": [
    "In a continuous model, we write the rate of change in the\n",
    "form of a derivative: \n",
    "\n",
    "$$\\frac{dx}{dt} = c$$ \n",
    "\n",
    "This way of representing the model is a *differential equation*, which is an equation that involves at least one derivative (see <http://modsimpy.com/diffeq>).\n",
    "To solve this equation, we multiply both sides by $dt$: \n",
    "\n",
    "$$dx = c dt$$ \n",
    "\n",
    "And then integrate both sides: \n",
    "\n",
    "$$x(t) = c t + x_0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "covered-defense",
   "metadata": {},
   "source": [
    "Similarly, we can write the proportional growth model like this:\n",
    "\n",
    "$$\\frac{\\Delta x}{\\Delta t} = \\alpha x$$ \n",
    "\n",
    "And as a differential equation\n",
    "like this: \n",
    "\n",
    "$$\\frac{dx}{dt} = \\alpha x$$ \n",
    "\n",
    "If we multiply both sides by\n",
    "$dt$ and divide by $x$, we get \n",
    "\n",
    "$$\\frac{1}{x}~dx = \\alpha~dt$$ \n",
    "\n",
    "Now we\n",
    "integrate both sides, yielding: \n",
    "\n",
    "$$\\ln x = \\alpha t + K$$ \n",
    "\n",
    "where $\\ln$ is the natural logarithm and $K$ is the constant of integration."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "underlying-visitor",
   "metadata": {},
   "source": [
    "Exponentiating both sides, we have\n",
    "\n",
    "$$\\exp(\\ln(x)) = \\exp(\\alpha t + K)$$ \n",
    "\n",
    "The exponential function can be written $\\exp(x)$ or $e^x$. In this book I use the first form because it resembles the Python code.\n",
    "We can rewrite the previous equation as\n",
    "\n",
    "$$x = \\exp(\\alpha t) \\exp(K)$$ \n",
    "\n",
    "Since $K$ is an arbitrary constant,\n",
    "$\\exp(K)$ is also an arbitrary constant, so we can write\n",
    "\n",
    "$$x = C \\exp(\\alpha t)$$ \n",
    "\n",
    "where $C = \\exp(K)$. There are many solutions to this differential equation, with different values of $C$. The particular solution we want is the one that has the value $x_0$ when $t=0$.\n",
    "\n",
    "When $t=0$, $x(t) = C$, so $C = x_0$ and the solution we want is\n",
    "\n",
    "$$x(t) = x_0 \\exp(\\alpha t)$$ \n",
    "\n",
    "If you would like to see this derivation done more carefully, you might like this video:\n",
    "<http://modsimpy.com/khan1>."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "czech-scratch",
   "metadata": {},
   "source": [
    "## Analysis and Simulation\n",
    "\n",
    "Once you have designed a model, there are generally two ways to proceed: simulation and analysis. Simulation often comes in the form of a computer program that models changes in a system over time, like births and deaths, or bikes moving from place to place. Analysis often comes in the form of algebra and calculus; that is, symbolic manipulation using mathematical notation.\n",
    "\n",
    "Analysis and simulation have different capabilities and limitations.\n",
    "Simulation is generally more versatile; it is easy to add and remove parts of a program and test many versions of a model, as we have done in the previous examples.\n",
    "\n",
    "But there are several things we can do with analysis that are harder or impossible with simulations:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "embedded-miniature",
   "metadata": {},
   "source": [
    "-   With analysis we can sometimes compute, exactly and efficiently, a\n",
    "    value that we could only approximate, less efficiently, with\n",
    "    simulation. For example, in the quadratic model we plotted net growth versus population and saw it crosses through zero when the population is\n",
    "    near 14 billion.  We could estimate the crossing point using a\n",
    "    numerical search algorithm (more about that later). But with a bit of algebra, we derived the general result that $K=-\\alpha/\\beta$.\n",
    "\n",
    "-   Analysis sometimes provides \"computational shortcuts\", that is, the\n",
    "    ability to jump forward in time to compute the state of a system\n",
    "    many time steps in the future without computing the intervening\n",
    "    states.\n",
    "\n",
    "-   We can use analysis to state and prove generalizations about models;\n",
    "    for example, we might prove that certain results will always or\n",
    "    never occur. With simulations, we can show examples and sometimes\n",
    "    find counterexamples, but it is hard to write proofs.\n",
    "\n",
    "-   Analysis can provide insight into models and the systems they\n",
    "    describe; for example, sometimes we can identify\n",
    "    qualitatively different ways the system can behave and key parameters that control\n",
    "    those behaviors."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "impressive-planning",
   "metadata": {},
   "source": [
    "When people see what analysis can do, they sometimes get drunk with\n",
    "power and imagine that it gives them a special ability to see past the veil of the material world and discern the laws of mathematics that govern the universe. When they analyze a model of a physical system, they talk about \"the math behind it\" as if our world is the mere shadow of a world of ideal mathematical entities (I am not making this up; see <http://modsimpy.com/plato>.).\n",
    "\n",
    "This is, of course, nonsense. Mathematical notation is a language\n",
    "designed by humans for a purpose, specifically to facilitate symbolic\n",
    "manipulations like algebra. Similarly, programming languages are\n",
    "designed for a purpose, specifically to represent computational ideas\n",
    "and run programs.\n",
    "\n",
    "Each of these languages is good for the purposes it was designed for and less good for other purposes. But they are often complementary, and one of the goals of this book is to show how they can be used together."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "graduate-conducting",
   "metadata": {},
   "source": [
    "## Analysis with WolframAlpha\n",
    "\n",
    "Until recently, most analysis was done by rubbing graphite on wood\n",
    "pulp, a process that is laborious and error-prone. A useful\n",
    "alternative is symbolic computation. If you have used services like\n",
    "WolframAlpha, you have used symbolic computation.\n",
    "\n",
    "For example, if you go to <https://www.wolframalpha.com/> and enter\n",
    "\n",
    "```\n",
    "df(t) / dt = alpha f(t)\n",
    "```\n",
    "\n",
    "WolframAlpha infers that `f(t)` is a function of `t` and `alpha` is a\n",
    "parameter; it classifies the query as a \"first-order linear ordinary\n",
    "differential equation\", and reports the general solution:\n",
    "\n",
    "$$f(t) = c_1 \\exp(\\alpha t)$$ \n",
    "\n",
    "If you add a second equation to specify the initial condition:\n",
    "\n",
    "```\n",
    "df(t) / dt = alpha f(t),  f(0) = p_0\n",
    "```\n",
    "\n",
    "WolframAlpha reports the particular solution:\n",
    "\n",
    "$$f(t) = p_0 \\exp(\\alpha t)$$\n",
    "\n",
    "WolframAlpha is based on Mathematica, a powerful programming language\n",
    "designed specifically for symbolic computation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "beginning-tamil",
   "metadata": {},
   "source": [
    "## Analysis with SymPy\n",
    "\n",
    "Python has a library called SymPy that provides symbolic computation\n",
    "tools similar to Mathematica. They are not as easy to use as\n",
    "WolframAlpha, but they have some other advantages.\n",
    "\n",
    "To use it, we'll define `Symbol` objects that represent names of variables and functions.\n",
    "The `symbols` function takes a string and returns `Symbol` objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "flush-balance",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import symbols\n",
    "\n",
    "t = symbols('t')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "faced-praise",
   "metadata": {},
   "source": [
    "Now when we use `t`, Python treats it like a variable name rather than a specific number. \n",
    "For example, if we use `t` as part of an expression, like this,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ancient-anderson",
   "metadata": {},
   "outputs": [],
   "source": [
    "expr = t + 1\n",
    "expr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "industrial-confidence",
   "metadata": {},
   "source": [
    "Python doesn't try to perform numerical addition; rather, it creates a\n",
    "new `Symbol` that represents the sum of `t` and `1`. We can evaluate the\n",
    "sum using `subs`, which substitutes a value for a symbol. This example\n",
    "substitutes 2 for `t`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "familiar-oakland",
   "metadata": {},
   "outputs": [],
   "source": [
    "expr.subs(t, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "progressive-celebration",
   "metadata": {},
   "source": [
    "Functions in SymPy are represented by a special kind of `Symbol`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "incomplete-sleep",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import Function\n",
    "\n",
    "f = Function('f')\n",
    "f"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "legislative-dispatch",
   "metadata": {},
   "source": [
    "Now if we write `f(t)`, we get an object that represents the evaluation of a function, $f$, at a value, $t$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "inside-sally",
   "metadata": {},
   "outputs": [],
   "source": [
    "f(t)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "interior-index",
   "metadata": {},
   "source": [
    "But again SymPy doesn't actually\n",
    "try to evaluate it."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "indian-trout",
   "metadata": {},
   "source": [
    "## Differential Equations In SymPy\n",
    "\n",
    "SymPy provides a function, `diff`, that can differentiate a function. We can apply it to `f(t)` like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "requested-program",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import diff\n",
    "\n",
    "dfdt = diff(f(t), t)\n",
    "dfdt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "treated-commissioner",
   "metadata": {},
   "source": [
    "The result is a `Symbol` that represents the derivative of `f` with\n",
    "respect to `t`. But again, SymPy doesn't try to compute the derivative\n",
    "yet.\n",
    "\n",
    "To represent a differential equation, we use `Eq`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "incomplete-parker",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import Eq\n",
    "\n",
    "alpha = symbols('alpha')\n",
    "eq1 = Eq(dfdt, alpha*f(t))\n",
    "eq1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "concerned-asthma",
   "metadata": {},
   "source": [
    "The result is an object that represents an equation.  Now\n",
    "we can use `dsolve` to solve this differential equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "described-overhead",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import dsolve\n",
    "\n",
    "solution_eq = dsolve(eq1)\n",
    "solution_eq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "inside-medicine",
   "metadata": {},
   "source": [
    "The result is the *general solution*, which still contains an unspecified constant, $C_1$.\n",
    "To get the *particular solution* where $f(0) = p_0$, we substitute $p_0$ for `C1`. \n",
    "First, we have to tell Python that `C1` is a symbol."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "material-voice",
   "metadata": {},
   "outputs": [],
   "source": [
    "C1 = symbols('C1')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "serious-license",
   "metadata": {},
   "source": [
    "Now we can substitute the value of $p_0$ for `C1`.\n",
    "For example, if $p_0$ is 1000:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "supposed-domestic",
   "metadata": {},
   "outputs": [],
   "source": [
    "particular = solution_eq.subs(C1, 1000)\n",
    "particular"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "modern-raleigh",
   "metadata": {},
   "source": [
    "When $t=0$, the value of $f(0)$ is $p_0$, which confirms that this is the solution we want."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "394e9410",
   "metadata": {},
   "outputs": [],
   "source": [
    "particular.subs(t, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "negative-chance",
   "metadata": {},
   "source": [
    "## Solving the Quadratic Growth Model\n",
    "\n",
    "To solve the quadratic growth curve, we'll use the `r, K` parameterization, so we'll need two more symbols:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "political-chinese",
   "metadata": {},
   "outputs": [],
   "source": [
    "r, K = symbols('r K')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "virtual-temperature",
   "metadata": {},
   "source": [
    "Now we can write the differential equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "lined-queen",
   "metadata": {},
   "outputs": [],
   "source": [
    "eq2 = Eq(diff(f(t), t), r * f(t) * (1 - f(t)/K))\n",
    "eq2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "relative-zoning",
   "metadata": {},
   "source": [
    "And solve it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "printable-typing",
   "metadata": {},
   "outputs": [],
   "source": [
    "solution_eq = dsolve(eq2)\n",
    "solution_eq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "convertible-simple",
   "metadata": {},
   "source": [
    "The result, `solution_eq`, contains `rhs`, which is the right-hand side of the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "final-treatment",
   "metadata": {},
   "outputs": [],
   "source": [
    "general = solution_eq.rhs\n",
    "general"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "phantom-melbourne",
   "metadata": {},
   "source": [
    "We can evaluate the right-hand side at $t=0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "freelance-admission",
   "metadata": {},
   "outputs": [],
   "source": [
    "at_0 = general.subs(t, 0)\n",
    "at_0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "precise-radar",
   "metadata": {},
   "source": [
    "Now we want to find the value of `C1` that makes `f(0) = p_0`.\n",
    "\n",
    "So we'll create the equation `at_0 = p_0` and solve for `C1`.  Because this is just an algebraic equation, not a differential equation, we use `solve`, not `dsolve`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "orange-glasgow",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import solve\n",
    "\n",
    "p_0 = symbols('p_0')\n",
    "solutions = solve(Eq(at_0, p_0), C1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aggressive-annual",
   "metadata": {},
   "source": [
    "The result from `solve` is a list of solutions.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "expensive-palace",
   "metadata": {},
   "outputs": [],
   "source": [
    "type(solutions), len(solutions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "periodic-bikini",
   "metadata": {},
   "source": [
    "In this case, there is only one solution, but we still get a list, so we have to use the bracket operator, `[0]`, to select the first one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "worthy-sleeve",
   "metadata": {},
   "outputs": [],
   "source": [
    "value_of_C1 = solutions[0]\n",
    "value_of_C1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "lightweight-pickup",
   "metadata": {},
   "source": [
    "Now in the general solution, we want to replace `C1` with the value of `C1` we just figured out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "aggressive-queue",
   "metadata": {},
   "outputs": [],
   "source": [
    "particular = general.subs(C1, value_of_C1)\n",
    "particular"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "electric-albania",
   "metadata": {},
   "source": [
    "The result is complicated, but SymPy provides a function that tries to simplify it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "funded-tolerance",
   "metadata": {},
   "outputs": [],
   "source": [
    "simpler = particular.simplify()\n",
    "simpler"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "scenic-preparation",
   "metadata": {},
   "source": [
    "This function is called the *logistic growth curve*; see\n",
    "<http://modsimpy.com/logistic>. In the context of growth models, the\n",
    "logistic function is often written like this:\n",
    "\n",
    "$$f(t) = \\frac{K}{1 + A \\exp(-rt)}$$ \n",
    "\n",
    "where $A = (K - p_0) / p_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "burning-sherman",
   "metadata": {
    "tags": []
   },
   "source": [
    "We can use SymPy to confirm that these two forms are equivalent.  First we represent the alternative version of the logistic function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "institutional-finish",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "A = (K - p_0) / p_0\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "noble-auditor",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sympy import exp\n",
    "\n",
    "logistic = K / (1 + A * exp(-r*t))\n",
    "logistic"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brave-conditions",
   "metadata": {
    "tags": []
   },
   "source": [
    "To see whether two expressions are equivalent, we can check whether their difference simplifies to 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "industrial-administrator",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "(particular - logistic).simplify()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "therapeutic-topic",
   "metadata": {
    "tags": []
   },
   "source": [
    "This test only works one way: if SymPy says the difference reduces to 0, the expressions are definitely equivalent (and not just numerically close).\n",
    "\n",
    "But if SymPy can't find a way to simplify the result to 0, that doesn't necessarily mean there isn't one.  Testing whether two expressions are equivalent is a surprisingly hard problem; in fact, there is no algorithm that can solve it in general."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "flying-bermuda",
   "metadata": {
    "tags": []
   },
   "source": [
    "If you use SymPy to compute an expression, and then want to evaluate that expression in Python, SymPy provides a function called `pycode` that generates Python code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "bibliographic-sarah",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sympy.printing.pycode import pycode\n",
    "\n",
    "pycode(simpler)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cutting-access",
   "metadata": {},
   "source": [
    "If you would like to see this differential equation solved by hand, you might like this video: <http://modsimpy.com/khan2>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "selective-calvin",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "In this chapter we wrote the growth models from the previous chapters in terms of difference and differential equations. We solved some of these equations by hand; for others, we used WolframAlpha and SymPy.\n",
    "\n",
    "What I called the \"constant growth\" model is more commonly called *linear growth* because the solution is a line.  If we model time as continuous, the solution is\n",
    "\n",
    "$$f(t) = p_0 + c t$$\n",
    "\n",
    "where $c$ is net annual growth.\n",
    "\n",
    "Similarly, the proportional growth model is usually called *exponential growth* because the solution is an exponential function:\n",
    "\n",
    "$$f(t) = p_0 \\exp{\\alpha t}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "outstanding-launch",
   "metadata": {},
   "source": [
    "Finally, the quadratic growth model is called *logistic growth* because the solution is a logistic function:\n",
    "\n",
    "$$f(t) = \\frac{K}{1 + A \\exp(-rt)}$$ \n",
    "\n",
    "where $A = (K - p_0) / p_0$.\n",
    "\n",
    "I avoided these terms until now because they are based on results we had not derived yet.\n",
    "\n",
    "With that, we are done modeling world population growth.\n",
    "The next chapter presents case studies where you can apply the tools we have learned so far."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "offshore-passage",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "This chapter is available as a Jupyter notebook where you can read the text, run the code, and work on the exercises. \n",
    "You can access the notebooks at <https://allendowney.github.io/ModSimPy/>."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ideal-ecology",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    " Use SymPy to solve the quadratic growth equation using the alternative parameterization\n",
    "\n",
    "$$ \\frac{df(t)}{dt} = \\alpha f(t) + \\beta f^2(t) $$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "internal-knock",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "smooth-sunglasses",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "fundamental-arlington",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "desirable-alpha",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "coupled-grounds",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "looking-ground",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "anticipated-paragraph",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "residential-spanking",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "\n",
    "  Use [WolframAlpha](https://www.wolframalpha.com/) to solve the quadratic growth model, using either or both forms of parameterization:\n",
    "\n",
    "    df(t) / dt = alpha f(t) + beta f(t)^2\n",
    "\n",
    "or\n",
    "\n",
    "    df(t) / dt = r f(t) (1 - f(t)/K)\n",
    "\n",
    "Find the general solution and also the particular solution where `f(0) = p_0`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "purple-traffic",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}