{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# App.E: Programming of differential equations\n", "\n", " **Hans Petter Langtangen**, Simula Research Laboratory and University of Oslo, Dept. of Informatics\n", "\n", "Date: **Aug 21, 2016**\n", "\n", "# How to solve any ordinary scalar differential equation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "{\\color{red}u'(t)} &= \\alpha {\\color{red}u(t)}(1 - R^{-1}{\\color{red}u(t)})\\\\\n", "{\\color{red}u(0)} &= U_0\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "
\n", "\n", "\n", "\n", "\n", "\n", "## Examples on scalar differential equations (ODEs)\n", "\n", "**Terminology:**\n", "\n", " * *Scalar ODE*: a single ODE, one unknown function\n", "\n", " * *Vector ODE* or *systems of ODEs*: several ODEs, several unknown functions\n", "\n", "\n", "\n", "**Examples:**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "u'&=\\alpha u\\quad\\hbox{exponential growth}\\\\\n", "u'&=\\alpha u\\left( 1-\\frac{u}{R}\\right)\\quad\\hbox{logistic growth}\\\\\n", "u' + b|u|u &= g\\quad\\hbox{falling body in fluid}\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## We shall write an ODE in a generic form: $u'=f(u,t)$\n", "\n", " * Our methods and software should be applicable to *any* ODE\n", "\n", " * Therefore we need an abstract notation for an arbitrary ODE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u^{\\prime}(t) = f(u(t), t)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The three ODEs on the last slide correspond to" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "f(u,t) &= \\alpha u,\\quad\\hbox{exponential growth}\\\\\n", "f(u,t) &= \\alpha u\\left( 1-\\frac{u}{R}\\right),\\quad\\hbox{logistic growth}\\\\\n", "f(u,t) &= -b|u|u + g,\\quad\\hbox{body in fluid}\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our task: write functions and classes that take $f$ as input and\n", "produce $u$ as output\n", "\n", "\n", "\n", "## What is the $f(u,t)$?\n", "\n", "**Problem:**\n", "\n", "Given an ODE," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\sqrt{u}u' - \\alpha(t) u^{3/2}(1-\\frac{u}{R(t)}) = 0,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "what is the $f(u,t)$?\n", "\n", "\n", "\n", "**Solution:**\n", "\n", "The target form is $u'=f(u,t)$, so we need to isolate $u'$ on the\n", "left-hand side:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u' = \\underbrace{\\alpha(t) u(1-\\frac{u}{R(t)})}_{f(u,t)}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Such abstract $f$ functions are widely used in mathematics\n", "\n", "**We can make generic software for:**\n", "\n", "\n", " * Numerical differentiation: $f'(x)$\n", "\n", " * Numerical integration: $\\int_a^b f(x)dx$\n", "\n", " * Numerical solution of algebraic equations: $f(x)=0$\n", "\n", "\n", "\n", "**Applications:**\n", "\n", "1. $\\frac{d}{dx} x^a\\sin (wx)$: $f(x)=x^a\\sin (wx)$\n", "\n", "2. $\\int_{-1}^1 (x^2\\tanh^{-1}x - (1+x^2)^{-1})dx$: $f(x)=x^2\\tanh^{-1}x - (1+x^2)^{-1}$, $a=-1$, $b=1$\n", "\n", "3. Solve $x^4\\sin x = \\tan x$: $f(x)=x^4\\sin x - \\tan x$\n", "\n", "\n", "\n", "## We use finite difference approximations to derivatives to turn an ODE into a difference equation\n", "\n", "**$u'=f(u,t)$.**\n", "\n", "Assume we have computed $u$ at discrete time points\n", "$t_0,t_1,\\ldots,t_k$. At $t_k$ we have the ODE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u'(t_k) = f(u(t_k),t_k)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Approximate $u'(t_k)$ by a forward finite difference," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u'(t_k)\\approx \\frac{u(t_{k+1})-u(t_k)}{\\Delta t}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Insert in the ODE at $t=t_k$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{u(t_{k+1})-u(t_k)}{\\Delta t} = f(u(t_k),t_k)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Terms with $u(t_k)$ are known, and this is an algebraic (difference) equation\n", "for $u(t_{k+1})$\n", "\n", "\n", "\n", "## The Forward Euler (or Euler's) method; idea\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## The Forward Euler (or Euler's) method; idea\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## The Forward Euler (or Euler's) method; mathematics\n", "\n", "Solving with respect to $u(t_{k+1})$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u(t_{k+1}) = u(t_k) + \\Delta t f(u(t_k), t_k)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is a very simple formula that we can use repeatedly for\n", "$u(t_1)$, $u(t_2)$, $u(t_3)$ and so forth.\n", "\n", "\n", "\n", "**Difference equation notation:**\n", "\n", "Let $u_k$ denote the numerical approximation to the exact solution $u(t)$\n", "at $t=t_k$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u_{k+1} = u_k + \\Delta t f(u_k,t_k)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is an ordinary difference equation we can solve!\n", "\n", "\n", "\n", "\n", "\n", "## Illustration of the forward finite difference\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Let's apply the method!\n", "\n", "**Problem: The world's simplest ODE.**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u' = u,\\quad t\\in (0,T]\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solve for $u$ at $t=t_k=k\\Delta t$, $k=0,1,2,\\ldots,t_n$,\n", "$t_0=0$, $t_n=T$\n", "\n", "\n", "\n", "**Forward Euler method:**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u_{k+1} = u_k + \\Delta t\\, f(u_k,t_k)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution by hand:**\n", "\n", "What is $f$? $f(u,t)=u$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u_{k+1} = u_k + \\Delta t f(u_k,t_k) = u_k + \\Delta t u_k = (1+\\Delta t)u_k\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First step:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u_1 = (1+\\Delta t) u_0\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "but what is $u_0$?\n", "\n", "\n", "\n", "## An ODE needs an initial condition: $u(0)=U_0$\n", "\n", "**Numerics:**\n", "\n", "Any ODE $u'=f(u,t)$ *must* have an initial condition $u(0)=U_0$, with\n", "known $U_0$, otherwise we cannot start the method!\n", "\n", "\n", "\n", "**Mathematics:**\n", "\n", "In mathematics: $u(0)=U_0$ must be specified to get a unique solution.\n", "\n", "\n", "\n", "**Example:**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u'=u\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solution: $u=Ce^t$ for any constant $C$. Say $u(0)=U_0$: $u=U_0e^t$.\n", "\n", "\n", "\n", "## We continue solution by hand\n", "\n", "Say $U_0=2$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "u_1 &= (1+\\Delta t) u_0 = (1+\\Delta t) U_0 = (1+\\Delta t)2 \\\\\n", "u_2 &= (1+\\Delta t) u_1 = (1+\\Delta t) (1+\\Delta t)2 = 2(1+\\Delta t)^2\\\\\n", "u_3 &= (1+\\Delta t) u_2 = (1+\\Delta t) 2(1+\\Delta t)^2 = 2(1+\\Delta t)^3\\\\\n", "u_4 &= (1+\\Delta t) u_3 = (1+\\Delta t) 2(1+\\Delta t)^3 = 2(1+\\Delta t)^4\\\\\n", "u_5 &= (1+\\Delta t) u_4 = (1+\\Delta t) 2(1+\\Delta t)^4 = 2(1+\\Delta t)^5\\\\\n", "\\vdots &= \\vdots\\\\\n", "u_k &= 2(1+\\Delta t)^k\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Actually, we found a formula for $u_k$! No need to program...\n", "\n", "\n", "\n", "## How accurate is our numerical method?\n", "\n", " * Exact solution: $u(t)=2e^t$, $u(t_k)=2e^{k\\Delta t} = 2(e^{\\Delta t})^k$\n", "\n", " * Numerical solution: $u_k = 2(1+\\Delta t)^k$\n", "\n", "When going from $t_k$ to $t_{k+1}$, the solution is amplified by a factor:\n", "\n", " * Exact: $u(t_{k+1}) = e^{\\Delta t}u(t_k)$\n", "\n", " * Numerical: $u_{k+1} = (1+\\Delta t)u_k$\n", "\n", "Using Taylor series for $e^x$ we see that" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "e^{\\Delta t} - (1+\\Delta t) = 1 + \\Delta t + \\frac{\\Delta t^2}{2}\n", "+ frac{\\Delta t^3}{6} + \\cdots - (1+\\Delta t) = frac{\\Delta t^3}{6}\n", "+ {\\cal O}(\\Delta t^4)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This error approaches 0 as $\\Delta t\\rightarrow 0$.\n", "\n", "\n", "\n", "## What about the general case $u'=f(u,t)$?\n", "\n", "Given any $U_0$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "u_1 &= u_0 + \\Delta t f(u_0, t_0)\\\\\n", "u_2 &= u_1 + \\Delta t f(u_1, t_1)\\\\\n", "u_3 &= u_2 + \\Delta t f(u_2, t_2)\\\\\n", "u_4 &= u_3 + \\Delta t f(u_3, t_3)\\\\\n", "&\\vdots\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "No general formula in this case...\n", "\n", "\n", "\n", "**Rule of thumb:**\n", "\n", "When hand calculations get boring, let's program!\n", "\n", "\n", "\n", "## We start with a specialized program for $u'=u$, $u(0)=U_0$\n", "\n", "**Algorithm:**\n", "\n", "Given $\\Delta t$ (`dt`) and $n$\n", "\n", " * Create arrays `t` and `u` of length $n+1$\n", "\n", " * Set initial condition: `u[0]` = $U_0$, `t[0]=0`\n", "\n", " * For $k=0,1,2,\\ldots,n-1$:\n", "\n", " * `t[k+1] = t[k] + dt`\n", "\n", " * `u[k+1] = (1 + dt)*u[k]`\n", "\n", "\n", "\n", "## We start with a specialized program for $u'=u$, $u(0)=U_0$\n", "\n", "**Program:**" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n", "import sys\n", "\n", "dt = float(sys.argv[1])\n", "U0 = 1\n", "T = 4\n", "n = int(T/dt)\n", "\n", "t = np.zeros(n+1)\n", "u = np.zeros(n+1)\n", "\n", "t[0] = 0\n", "u[0] = U0\n", "for k in range(n):\n", " t[k+1] = t[k] + dt\n", " u[k+1] = (1 + dt)*u[k]\n", "\n", "# plot u against t" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The solution if we plot $u$ against $t$\n", "\n", "$\\Delta t = 0.4$ and $\\Delta t=0.2$:\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## The algorithm for the general ODE $u'=f(u,t)$\n", "\n", "**Algorithm:**\n", "\n", "Given $\\Delta t$ (`dt`) and $n$\n", "\n", " * Create arrays `t` and `u` of length $n+1$\n", "\n", " * Create array `u` to hold $u_k$ and\n", "\n", " * Set initial condition: `u[0]` = $U_0$, `t[0]=0`\n", "\n", " * For $k=0,1,2,\\ldots,n-1$:\n", "\n", " * `u[k+1] = u[k] + dt*f(u[k], t[k])` (the only change!)\n", "\n", " * `t[k+1] = t[k] + dt`\n", "\n", "\n", "\n", "## Implementation of the general algorithm for $u'=f(u,t)$\n", "\n", "**General function:**" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def ForwardEuler(f, U0, T, n):\n", " \"\"\"Solve u'=f(u,t), u(0)=U0, with n steps until t=T.\"\"\"\n", " import numpy as np\n", " t = np.zeros(n+1)\n", " u = np.zeros(n+1) # u[k] is the solution at time t[k]\n", "\n", " u[0] = U0\n", " t[0] = 0\n", " dt = T/float(n)\n", "\n", " for k in range(n):\n", " t[k+1] = t[k] + dt\n", " u[k+1] = u[k] + dt*f(u[k], t[k])\n", "\n", " return u, t" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Magic:**\n", "\n", "This simple function can solve any ODE (!)\n", "\n", "\n", "\n", "## Example on using the function\n", "\n", "**Mathematical problem:**\n", "\n", "Solve $u'=u$, $u(0)=1$, for $t\\in [0,4]$,\n", "with $\\Delta t = 0.4$ \n", "Exact solution: $u(t)=e^t$.\n", "\n", "\n", "\n", "**Basic code:**" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def f(u, t):\n", " return u\n", "\n", "U0 = 1\n", "T = 3\n", "n = 30\n", "u, t = ForwardEuler(f, U0, T, n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Compare exact and numerical solution:**" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "from scitools.std import plot, exp\n", "u_exact = exp(t)\n", "plot(t, u, 'r-', t, u_exact, 'b-',\n", " xlabel='t', ylabel='u', legend=('numerical', 'exact'),\n", " title=\"Solution of the ODE u'=u, u(0)=1\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Now you can solve any ODE!\n", "\n", "**Recipe:**\n", "\n", " * Identify $f(u,t)$ in your ODE\n", "\n", " * Make sure you have an initial condition $U_0$\n", "\n", " * Implement the $f(u,t)$ formula in a Python function `f(u, t)`\n", "\n", " * Choose $\\Delta t$ or no of steps $n$\n", "\n", " * Call `u, t = ForwardEuler(f, U0, T, n)`\n", "\n", " * `plot(t, u)`\n", "\n", "\n", "\n", "**Warning:**\n", "\n", "The Forward Euler method may give very inaccurate solutions if\n", "$\\Delta t$ is not sufficiently small. For some problems\n", "(like $u''+u=0$) other methods should be used.\n", "\n", "\n", "\n", "\n", "## Let us make a class instead of a function for solving ODEs\n", "\n", "**Usage of the class:**" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "method = ForwardEuler(f, dt)\n", "method.set_initial_condition(U0, t0)\n", "u, t = method.solve(T)\n", "plot(t, u)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**How?**\n", "\n", " * Store $f$, $\\Delta t$, and the sequences $u_k$, $t_k$ as attributes\n", "\n", " * Split the steps in the `ForwardEuler` function into four methods:\n", "\n", " * the constructor (`__init__`)\n", "\n", " * `set_initial_condition` for $u(0)=U_0$\n", "\n", " * `solve` for running the numerical time stepping\n", "\n", " * `advance` for isolating the numerical updating formula \n", " (new numerical methods just need a different `advance` method,\n", " the rest is the same)\n", "\n", "\n", "\n", "## The code for a class for solving ODEs (part 1)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "\n", "import numpy as np\n", "\n", "class ForwardEuler_v1:\n", " def __init__(self, f, dt):\n", " self.f, self.dt = f, dt\n", "\n", " def set_initial_condition(self, U0):\n", " self.U0 = float(U0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The code for a class for solving ODEs (part 2)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [], "source": [ "\n", "class ForwardEuler_v1:\n", " ...\n", " def solve(self, T):\n", " \"\"\"Compute solution for 0 <= t <= T.\"\"\"\n", " n = int(round(T/self.dt)) # no of intervals\n", " self.u = np.zeros(n+1)\n", " self.t = np.zeros(n+1)\n", " self.u[0] = float(self.U0)\n", " self.t[0] = float(0)\n", "\n", " for k in range(self.n):\n", " self.k = k\n", " self.t[k+1] = self.t[k] + self.dt\n", " self.u[k+1] = self.advance()\n", " return self.u, self.t\n", "\n", " def advance(self):\n", " \"\"\"Advance the solution one time step.\"\"\"\n", " # Create local variables to get rid of \"self.\" in\n", " # the numerical formula\n", " u, dt, f, k, t = self.u, self.dt, self.f, self.k, self.t\n", "\n", " unew = u[k] + dt*f(u[k], t[k])\n", " return unew" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Alternative class code for solving ODEs (part 1)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [], "source": [ "\n", "# Idea: drop dt in the constructor.\n", "# Let the user provide all time points (in solve).\n", "\n", "class ForwardEuler:\n", " def __init__(self, f):\n", "\n", " # test that f is a function\n", " if not callable(f):\n", " raise TypeError('f is %s, not a function' % type(f))\n", " self.f = f\n", "\n", " def set_initial_condition(self, U0):\n", " self.U0 = float(U0)\n", "\n", " def solve(self, time_points):\n", " ..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Alternative class code for solving ODEs (part 2)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [], "source": [ "\n", "class ForwardEuler:\n", " ...\n", " def solve(self, time_points):\n", " \"\"\"Compute u for t values in time_points list.\"\"\"\n", " self.t = np.asarray(time_points)\n", " self.u = np.zeros(len(time_points))\n", "\n", " self.u[0] = self.U0\n", "\n", " for k in range(len(self.t)-1):\n", " self.k = k\n", " self.u[k+1] = self.advance()\n", " return self.u, self.t\n", "\n", " def advance(self):\n", " \"\"\"Advance the solution one time step.\"\"\"\n", " u, f, k, t = self.u, self.f, self.k, self.t\n", "\n", " dt = t[k+1] - t[k]\n", " unew = u[k] + dt*f(u[k], t[k])\n", " return unew" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Verifying the class implementation; mathematics\n", "\n", "**Mathematical problem:**\n", "\n", "Important result: the numerical method (and most others) will\n", "exactly reproduce $u$ if it is linear in $t$ (!):" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "u(t) &= at+b = 0.2t + 3\\\\\n", "h(t) &= u(t)\\\\\n", "u'(t) &=0.2 + (u-h(t))^4,\\quad u(0)=3,\\quad t\\in [0,3]\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This $u$ should be reproduced to machine precision for any $\\Delta t$.\n", "\n", "\n", "\n", "## Verifying the class implementation; implementation\n", "\n", "**Code:**" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def test_ForwardEuler_against_linear_solution():\n", " def f(u, t):\n", " return 0.2 + (u - h(t))**4\n", "\n", " def h(t):\n", " return 0.2*t + 3\n", "\n", " solver = ForwardEuler(f)\n", " solver.set_initial_condition(U0=3)\n", " dt = 0.4; T = 3; n = int(round(float(T)/dt))\n", " time_points = np.linspace(0, T, n+1)\n", " u, t = solver.solve(time_points)\n", " u_exact = h(t)\n", " diff = np.abs(u_exact - u).max()\n", " tol = 1E-14\n", " success = diff < tol\n", " assert success" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Using a class to hold the right-hand side $f(u,t)$\n", "\n", "**Mathematical problem:**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u^{\\prime}(t)=\\alpha u(t)\\left( 1-\\frac{u(t)}{R}\\right),\\quad u(0)=U_0,\\quad t\\in [0,40]\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Class for right-hand side $f(u,t)$:**" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [], "source": [ "\n", "class Logistic:\n", " def __init__(self, alpha, R, U0):\n", " self.alpha, self.R, self.U0 = alpha, float(R), U0\n", "\n", " def __call__(self, u, t): # f(u,t)\n", " return self.alpha*u*(1 - u/self.R)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Main program:**" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [], "source": [ "problem = Logistic(0.2, 1, 0.1)\n", "time_points = np.linspace(0, 40, 401)\n", "method = ForwardEuler(problem)\n", "method.set_initial_condition(problem.U0)\n", "u, t = method.solve(time_points)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Figure of the solution\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Numerical methods for ordinary differential equations\n", "\n", "**Forward Euler method:**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u_{k+1} = u_k + \\Delta t\\, f(u_k, t_k)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**4th-order Runge-Kutta method:**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u_{k+1} = u_k + {1\\over 6}\\left( K_1 + 2K_2 + 2K_3 + K_4\\right)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "K_1 &= \\Delta t\\,f(u_k, t_k)\\\\\n", "K_2 &= \\Delta t\\,f(u_k + {1\\over2}K_1, t_k + {1\\over2}\\Delta t)\\\\\n", "K_3 &= \\Delta t\\,f(u_k + {1\\over2}K_2, t_k + {1\\over2}\\Delta t)\\\\\n", "K_4 &= \\Delta t\\,f(u_k + K3, t_k + \\Delta t)\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And lots of other methods! How to program a wide collection of methods?\n", "Use object-oriented programming!\n", "\n", "\n", "\n", "## A superclass for ODE methods\n", "\n", "**Common tasks for ODE solvers:**\n", "\n", " * Store the solution $u_k$ and the corresponding time levels $t_k$, $k=0,1,2,\\ldots,n$\n", "\n", " * Store the right-hand side function $f(u,t)$\n", "\n", " * Set and store the initial condition\n", "\n", " * Run the loop over all time steps\n", "\n", "\n", "\n", "**Principles:**\n", "\n", " * Common data and functionality are placed in superclass `ODESolver`\n", "\n", " * Isolate the numerical updating formula in a method `advance`\n", "\n", " * Subclasses, e.g., `ForwardEuler`, just implement the specific numerical formula in `advance`\n", "\n", "\n", "\n", "## The superclass code" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [], "source": [ "class ODESolver:\n", " def __init__(self, f):\n", " self.f = f\n", "\n", " def advance(self):\n", " \"\"\"Advance solution one time step.\"\"\"\n", " raise NotImplementedError # implement in subclass\n", "\n", " def set_initial_condition(self, U0):\n", " self.U0 = float(U0)\n", "\n", " def solve(self, time_points):\n", " self.t = np.asarray(time_points)\n", " self.u = np.zeros(len(self.t))\n", " # Assume that self.t[0] corresponds to self.U0\n", " self.u[0] = self.U0\n", "\n", " # Time loop\n", " for k in range(n-1):\n", " self.k = k\n", " self.u[k+1] = self.advance()\n", " return self.u, self.t\n", "\n", " def advance(self):\n", " raise NotImplemtedError # to be impl. in subclasses" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Implementation of the Forward Euler method\n", "\n", "**Subclass code:**" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [], "source": [ "class ForwardEuler(ODESolver):\n", " def advance(self):\n", " u, f, k, t = self.u, self.f, self.k, self.t\n", "\n", " dt = t[k+1] - t[k]\n", " unew = u[k] + dt*f(u[k], t)\n", " return unew" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Application code for $u'-u=0$, $u(0)=1$, $t\\in [0,3]$, $\\Delta t=0.1$:**" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from ODESolver import ForwardEuler\n", "def test1(u, t):\n", " return u\n", "\n", "method = ForwardEuler(test1)\n", "method.set_initial_condition(U0=1)\n", "u, t = method.solve(time_points=np.linspace(0, 3, 31))\n", "plot(t, u)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The implementation of a Runge-Kutta method\n", "\n", "**Subclass code:**" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [], "source": [ "class RungeKutta4(ODESolver):\n", " def advance(self):\n", " u, f, k, t = self.u, self.f, self.k, self.t\n", "\n", " dt = t[k+1] - t[k]\n", " dt2 = dt/2.0\n", " K1 = dt*f(u[k], t)\n", " K2 = dt*f(u[k] + 0.5*K1, t + dt2)\n", " K3 = dt*f(u[k] + 0.5*K2, t + dt2)\n", " K4 = dt*f(u[k] + K3, t + dt)\n", " unew = u[k] + (1/6.0)*(K1 + 2*K2 + 2*K3 + K4)\n", " return unew" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Application code (same as for ForwardEuler):**" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from ODESolver import RungeKutta4\n", "def test1(u, t):\n", " return u\n", "\n", "method = RungeKutta4(test1)\n", "method.set_initial_condition(U0=1)\n", "u, t = method.solve(time_points=np.linspace(0, 3, 31))\n", "plot(t, u)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The user should be able to check intermediate solutions and terminate the time stepping\n", "\n", " * Sometimes a property of the solution determines when to stop the solution process: e.g., when $u < 10^{-7}\\approx 0$.\n", "\n", " * Extension: `solve(time_points, terminate)`\n", "\n", " * `terminate(u, t, step_no)` is called at every time step, is user-defined,\n", " and returns `True` when the time stepping should be terminated\n", "\n", " * Last computed solution is `u[step_no]` at time `t[step_no]`" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def terminate(u, t, step_no):\n", " eps = 1.0E-6 # small number\n", " return abs(u[step_no,0]) < eps # close enough to zero?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Systems of differential equations (vector ODE)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "u' &= v\\\\\n", "v' &= -u\\\\\n", "u(0) &= 1\\\\\n", "v(0) &= 0\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Example on a system of ODEs (vector ODE)\n", "\n", "Two ODEs with two unknowns $u(t)$ and $v(t)$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", " u'(t) &= v(t)\\\\\n", " v'(t) &= -u(t)\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Each unknown must have an initial condition, say" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u(0)=0,\\quad v(0)=1\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this case, one can derive the exact solution to be" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u(t)=\\sin (t),\\quad v(t)=\\cos (t)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Systems of ODEs appear frequently in physics, biology, finance, ...\n", "\n", "\n", "\n", "## The ODE system that is the final project in the course\n", "\n", "Model for spreading of a disease in a population:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "S' &= - \\beta SI\\\\\n", "I' &= \\beta SI -\\nu R\\\\\n", "R' &= \\nu I\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Initial conditions:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "S(0) &= S_0\\\\\n", "I(0) &= I_0\\\\\n", "R(0) &= 0\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Another example on a system of ODEs (vector ODE)\n", "\n", "Second-order ordinary differential equation, for a spring-mass system\n", "(from Newton's second law):" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "mu'' + \\beta u' + ku = 0,\\quad u(0)=U_0,\\ u'(0)=0\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can rewrite this as a system of two *first-order* equations, by\n", "introducing two new unknowns" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u^{(0)}(t) \\equiv u(t),\\quad u^{(1)}(t) \\equiv u'(t)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first-order system is then" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "{d\\over dt} u^{(0)}(t) &= u^{(1)}(t)\\\\\n", "{d\\over dt} u^{(1)}(t) &= -m^{-1}\\beta u^{(1)} - m^{-1}ku^{(0)}\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Initial conditions: $u^{(0)}(0) = U_0$, $u^{(1)}(0)=0$\n", "\n", "\n", "\n", "## Making a flexible toolbox for solving ODEs\n", "\n", " * For scalar ODEs we could make one general class hierarchy to solve\n", " \"all\" problems with a range of methods\n", "\n", " * Can we easily extend class hierarchy to systems of ODEs?\n", "\n", " * Yes!\n", "\n", " * The example here can easily be extended to professional code\n", " ([Odespy](https://github.com/hplgit/odespy))\n", "\n", "\n", "\n", "## Vector notation for systems of ODEs: unknowns and equations\n", "\n", "General software for any vector/scalar ODE demands a general\n", "mathematical notation. We introduce $n$ unknowns" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u^{(0)}(t), u^{(1)}(t), \\ldots, u^{(n-1)}(t)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "in a system of $n$ ODEs:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "{d\\over dt}u^{(0)} &= f^{(0)}(u^{(0)}, u^{(1)}, \\ldots, u^{(n-1)}, t)\\\\\n", "{d\\over dt}u^{(1)} &= f^{(1)}(u^{(0)}, u^{(1)}, \\ldots, u^{(n-1)}, t)\\\\\n", "\\vdots &=& \\vdots\\\\\n", "{d\\over dt}u^{(n-1)} &= f^{(n-1)}(u^{(0)}, u^{(1)}, \\ldots, u^{(n-1)}, t)\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Vector notation for systems of ODEs: vectors\n", "\n", "We can collect the $u^{(i)}(t)$ functions and right-hand side functions $f^{(i)}$ in vectors:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u = (u^{(0)}, u^{(1)}, \\ldots, u^{(n-1)})\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "f = (f^{(0)}, f^{(1)}, \\ldots, f^{(n-1)})\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first-order system can then be written" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u' = f(u, t),\\quad u(0) = U_0\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $u$ and $f$ are vectors and $U_0$ is a vector of initial conditions\n", "\n", "\n", "\n", "**The magic of this notation:**\n", "\n", "Observe that the notation makes a scalar ODE and a system look the same,\n", "and we can easily make Python code that can handle both cases within\n", "the same lines of code (!)\n", "\n", "\n", "\n", "\n", "## How to make class ODESolver work for systems of ODEs\n", "\n", " * Recall: ODESolver was written for a scalar ODE\n", "\n", " * Now we want it to work for a system $u'=f$, $u(0)=U_0$, where $u$, $f$ and $U_0$ are vectors (arrays)\n", "\n", " * What are the problems?\n", "\n", "Forward Euler applied to a system:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\underbrace{u_{k+1}}_{\\mbox{vector}} =\n", "\\underbrace{u_k}_{\\mbox{vector}} +\n", "\\Delta t\\, \\underbrace{f(u_k, t_k)}_{\\mbox{vector}}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In Python code:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [], "source": [ "unew = u[k] + dt*f(u[k], t)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where\n", "\n", " * `u` is a two-dim. array (`u[k]` is a row)\n", "\n", " * `f` is a function returning an array\n", " (all the right-hand sides $f^{(0)},\\ldots,f^{(n-1)}$)\n", "\n", "\n", "\n", "## The adjusted superclass code (part 1)\n", "\n", "**To make ODESolver work for systems:**\n", "\n", " * Ensure that `f(u,t)` returns an array. \n", " This can be done be a general adjustment in the superclass!\n", "\n", " * Inspect $U_0$ to see if it is a number or list/tuple and\n", " make corresponding `u` 1-dim or 2-dim array" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [], "source": [ "class ODESolver:\n", " def __init__(self, f):\n", " # Wrap user's f in a new function that always\n", " # converts list/tuple to array (or let array be array)\n", " self.f = lambda u, t: np.asarray(f(u, t), float)\n", "\n", " def set_initial_condition(self, U0):\n", " if isinstance(U0, (float,int)): # scalar ODE\n", " self.neq = 1 # no of equations\n", " U0 = float(U0)\n", " else: # system of ODEs\n", " U0 = np.asarray(U0)\n", " self.neq = U0.size # no of equations\n", " self.U0 = U0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The superclass code (part 2)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [], "source": [ "class ODESolver:\n", " ...\n", " def solve(self, time_points, terminate=None):\n", " if terminate is None:\n", " terminate = lambda u, t, step_no: False\n", "\n", " self.t = np.asarray(time_points)\n", " n = self.t.size\n", " if self.neq == 1: # scalar ODEs\n", " self.u = np.zeros(n)\n", " else: # systems of ODEs\n", " self.u = np.zeros((n,self.neq))\n", "\n", " # Assume that self.t[0] corresponds to self.U0\n", " self.u[0] = self.U0\n", "\n", " # Time loop\n", " for k in range(n-1):\n", " self.k = k\n", " self.u[k+1] = self.advance()\n", " if terminate(self.u, self.t, self.k+1):\n", " break # terminate loop over k\n", " return self.u[:k+2], self.t[:k+2]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "All subclasses from the scalar ODE works for systems as well\n", "\n", "\n", "\n", "## Example on how to use the general class hierarchy\n", "\n", "**Spring-mass system formulated as a system of ODEs:**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "mu'' + \\beta u' + ku = 0,\\quad u(0),\\ u'(0)\\hbox{ known}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "u^{(0)} &= u,\\quad u^{(1)} = u'\\\\\n", "u(t) &= (u^{(0)}(t), u^{(1)}(t))\\\\\n", "f(u,t)&=(u^{(1)}(t), -m^{-1}\\beta u^{(1)} - m^{-1}ku^{(0)})\\\\\n", "u'(t) &= f(u,t)\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Code defining the right-hand side:**" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def myf(u, t):\n", " # u is array with two components u[0] and u[1]:\n", " return [u[1],\n", " -beta*u[1]/m - k*u[0]/m]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Alternative implementation of the $f$ function via a class\n", "\n", "**Better (no global variables):**" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [], "source": [ "class MyF:\n", " def __init__(self, m, k, beta):\n", " self.m, self.k, self.beta = m, k, beta\n", "\n", " def __call__(self, u, t):\n", " m, k, beta = self.m, self.k, self.beta\n", " return [u[1], -beta*u[1]/m - k*u[0]/m]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Main program:**" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from ODESolver import ForwardEuler\n", "# initial condition:\n", "U0 = [1.0, 0]\n", "f = MyF(1.0, 1.0, 0.0) # u'' + u = 0 => u(t)=cos(t)\n", "solver = ForwardEuler(f)\n", "solver.set_initial_condition(U0)\n", "\n", "T = 4*pi; dt = pi/20; n = int(round(T/dt))\n", "time_points = np.linspace(0, T, n+1)\n", "u, t = solver.solve(time_points)\n", "\n", "# u is an array of [u0,u1] arrays, plot all u0 values:\n", "u0_values = u[:,0]\n", "u0_exact = cos(t)\n", "plot(t, u0_values, 'r-', t, u0_exact, 'b-')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Throwing a ball; ODE model\n", "\n", "Newton's 2nd law for a ball's trajectory through air leads to" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "{dx\\over dt} &= v_x\\\\\n", "{dv_x\\over dt} &= 0\\\\\n", "{dy\\over dt} &= v_y\\\\\n", "{dv_y\\over dt} &= -g\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Air resistance is neglected but can easily be added!\n", "\n", " * 4 ODEs with 4 unknowns:\n", "\n", " * the ball's position $x(t)$, $y(t)$\n", "\n", " * the velocity $v_x(t)$, $v_y(t)$\n", "\n", "\n", "\n", "## Throwing a ball; code\n", "\n", "**Define the right-hand side:**" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def f(u, t):\n", " x, vx, y, vy = u\n", " g = 9.81\n", " return [vx, 0, vy, -g]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Main program:**" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from ODESolver import ForwardEuler\n", "# t=0: prescribe x, y, vx, vy\n", "x = y = 0 # start at the origin\n", "v0 = 5; theta = 80*pi/180 # velocity magnitude and angle\n", "vx = v0*cos(theta)\n", "vy = v0*sin(theta)\n", "# Initial condition:\n", "U0 = [x, vx, y, vy]\n", "\n", "solver= ForwardEuler(f)\n", "solver.set_initial_condition(u0)\n", "time_points = np.linspace(0, 1.2, 101)\n", "u, t = solver.solve(time_points)\n", "\n", "# u is an array of [x,vx,y,vy] arrays, plot y vs x:\n", "x = u[:,0]; y = u[:,2]\n", "plot(x, y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Throwing a ball; results\n", "\n", "\n", "Comparison of exact and Forward Euler solutions\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Logistic growth model; ODE and code overview\n", "\n", "**Model:**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "u' = \\alpha u(1 - u/{\\color{red}R(t)}),\\quad u(0)=U_0\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$R$ is the maximum population size,\n", "which can vary with changes in the environment over time\n", "\n", "\n", "\n", "**Implementation features:**\n", "\n", " * Class Problem holds \"all physics\": $\\alpha$, $R(t)$, $U_0$, $T$ (end time), $f(u,t)$ in ODE\n", "\n", " * Class Solver holds \"all numerics\": $\\Delta t$, solution method; solves the problem and plots the solution\n", "\n", " * Solve for $t\\in [0,T]$ but terminate when $|u-R| < \\hbox{tol}$\n", "\n", "\n", "\n", "## Logistic growth model; class Problem ($f$)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [], "source": [ "class Problem:\n", " def __init__(self, alpha, R, U0, T):\n", " self.alpha, self.R, self.U0, self.T = alpha, R, U0, T\n", "\n", " def __call__(self, u, t):\n", " \"\"\"Return f(u, t).\"\"\"\n", " return self.alpha*u*(1 - u/self.R(t))\n", "\n", " def terminate(self, u, t, step_no):\n", " \"\"\"Terminate when u is close to R.\"\"\"\n", " tol = self.R*0.01\n", " return abs(u[step_no] - self.R) < tol\n", "\n", "problem = Problem(alpha=0.1, R=500, U0=2, T=130)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Logistic growth model; class Solver" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [], "source": [ "class Solver:\n", " def __init__(self, problem, dt,\n", " method=ODESolver.ForwardEuler):\n", " self.problem, self.dt = problem, dt\n", " self.method = method\n", "\n", " def solve(self):\n", " solver = self.method(self.problem)\n", " solver.set_initial_condition(self.problem.U0)\n", " n = int(round(self.problem.T/self.dt))\n", " t_points = np.linspace(0, self.problem.T, n+1)\n", " self.u, self.t = solver.solve(t_points,\n", " self.problem.terminate)\n", "\n", " def plot(self):\n", " plot(self.t, self.u)\n", "\n", "problem = Problem(alpha=0.1, U0=2, T=130,\n", " R=lambda t: 500 if t < 60 else 100)\n", "solver = Solver(problem, dt=1.)\n", "solver.solve()\n", "solver.plot()\n", "print 'max u:', solver.u.max()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Logistic growth model; results\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "" ] } ], "metadata": {}, "nbformat": 4, "nbformat_minor": 0 }