App.E: Programming of differential equations

Hans Petter Langtangen [1, 2]

[1] Simula Research Laboratory
[2] University of Oslo, Dept. of Informatics

Aug 21, 2016










How to solve any ordinary scalar differential equation

$$ \begin{align*} {\color{red}u'(t)} &= \alpha {\color{red}u(t)}(1 - R^{-1}{\color{red}u(t)})\\ {\color{red}u(0)} &= U_0 \end{align*} $$













Examples on scalar differential equations (ODEs)

Terminology:

Examples:

$$ \begin{align*} u'&=\alpha u\quad\hbox{exponential growth}\\ u'&=\alpha u\left( 1-\frac{u}{R}\right)\quad\hbox{logistic growth}\\ u' + b|u|u &= g\quad\hbox{falling body in fluid} \end{align*} $$









We shall write an ODE in a generic form: \( u'=f(u,t) \)

$$ u^{\prime}(t) = f(u(t), t)$$

The three ODEs on the last slide correspond to $$ \begin{align*} f(u,t) &= \alpha u,\quad\hbox{exponential growth}\\ f(u,t) &= \alpha u\left( 1-\frac{u}{R}\right),\quad\hbox{logistic growth}\\ f(u,t) &= -b|u|u + g,\quad\hbox{body in fluid} \end{align*} $$

Our task: write functions and classes that take \( f \) as input and produce \( u \) as output









What is the \( f(u,t) \)?

Problem:

Given an ODE, $$ \sqrt{u}u' - \alpha(t) u^{3/2}(1-\frac{u}{R(t)}) = 0,$$ what is the \( f(u,t) \)?

Solution:

The target form is \( u'=f(u,t) \), so we need to isolate \( u' \) on the left-hand side: $$ u' = \underbrace{\alpha(t) u(1-\frac{u}{R(t)})}_{f(u,t)} $$









Such abstract \( f \) functions are widely used in mathematics

We can make generic software for:

Applications:

  1. \( \frac{d}{dx} x^a\sin (wx) \): \( f(x)=x^a\sin (wx) \)
  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 \)
  3. Solve \( x^4\sin x = \tan x \): \( f(x)=x^4\sin x - \tan x \)









We use finite difference approximations to derivatives to turn an ODE into a difference equation

\( u'=f(u,t) \).

Assume we have computed \( u \) at discrete time points \( t_0,t_1,\ldots,t_k \). At \( t_k \) we have the ODE $$ u'(t_k) = f(u(t_k),t_k) $$

Approximate \( u'(t_k) \) by a forward finite difference, $$ u'(t_k)\approx \frac{u(t_{k+1})-u(t_k)}{\Delta t}$$

Insert in the ODE at \( t=t_k \): $$ \frac{u(t_{k+1})-u(t_k)}{\Delta t} = f(u(t_k),t_k) $$

Terms with \( u(t_k) \) are known, and this is an algebraic (difference) equation for \( u(t_{k+1}) \)









The Forward Euler (or Euler's) method; idea













The Forward Euler (or Euler's) method; idea













The Forward Euler (or Euler's) method; mathematics

Solving with respect to \( u(t_{k+1}) \) $$ u(t_{k+1}) = u(t_k) + \Delta t f(u(t_k), t_k)$$

This is a very simple formula that we can use repeatedly for \( u(t_1) \), \( u(t_2) \), \( u(t_3) \) and so forth.

Difference equation notation:

Let \( u_k \) denote the numerical approximation to the exact solution \( u(t) \) at \( t=t_k \). $$ u_{k+1} = u_k + \Delta t f(u_k,t_k)$$

This is an ordinary difference equation we can solve!









Illustration of the forward finite difference













Let's apply the method!

Problem: The world's simplest ODE.

$$ u' = u,\quad t\in (0,T] $$

Solve for \( u \) at \( t=t_k=k\Delta t \), \( k=0,1,2,\ldots,t_n \), \( t_0=0 \), \( t_n=T \)

Forward Euler method:

$$ u_{k+1} = u_k + \Delta t\, f(u_k,t_k)$$

Solution by hand:

What is \( f \)? \( f(u,t)=u \) $$ u_{k+1} = u_k + \Delta t f(u_k,t_k) = u_k + \Delta t u_k = (1+\Delta t)u_k$$

First step: $$ u_1 = (1+\Delta t) u_0$$ but what is \( u_0 \)?









An ODE needs an initial condition: \( u(0)=U_0 \)

Numerics:

Any ODE \( u'=f(u,t) \) must have an initial condition \( u(0)=U_0 \), with known \( U_0 \), otherwise we cannot start the method!

Mathematics:

In mathematics: \( u(0)=U_0 \) must be specified to get a unique solution.

Example:

$$ u'=u $$ Solution: \( u=Ce^t \) for any constant \( C \). Say \( u(0)=U_0 \): \( u=U_0e^t \).









We continue solution by hand

Say \( U_0=2 \): $$ \begin{align*} u_1 &= (1+\Delta t) u_0 = (1+\Delta t) U_0 = (1+\Delta t)2 \\ u_2 &= (1+\Delta t) u_1 = (1+\Delta t) (1+\Delta t)2 = 2(1+\Delta t)^2\\ u_3 &= (1+\Delta t) u_2 = (1+\Delta t) 2(1+\Delta t)^2 = 2(1+\Delta t)^3\\ u_4 &= (1+\Delta t) u_3 = (1+\Delta t) 2(1+\Delta t)^3 = 2(1+\Delta t)^4\\ u_5 &= (1+\Delta t) u_4 = (1+\Delta t) 2(1+\Delta t)^4 = 2(1+\Delta t)^5\\ \vdots &= \vdots\\ u_k &= 2(1+\Delta t)^k \end{align*} $$ Actually, we found a formula for \( u_k \)! No need to program...









How accurate is our numerical method?

When going from \( t_k \) to \( t_{k+1} \), the solution is amplified by a factor: Using Taylor series for \( e^x \) we see that $$ e^{\Delta t} - (1+\Delta t) = 1 + \Delta t + \frac{\Delta t^2}{2} + frac{\Delta t^3}{6} + \cdots - (1+\Delta t) = frac{\Delta t^3}{6} + {\cal O}(\Delta t^4)$$ This error approaches 0 as \( \Delta t\rightarrow 0 \).









What about the general case \( u'=f(u,t) \)?

Given any \( U_0 \): $$ \begin{align*} u_1 &= u_0 + \Delta t f(u_0, t_0)\\ u_2 &= u_1 + \Delta t f(u_1, t_1)\\ u_3 &= u_2 + \Delta t f(u_2, t_2)\\ u_4 &= u_3 + \Delta t f(u_3, t_3)\\ &\vdots \end{align*} $$ No general formula in this case...

Rule of thumb:

When hand calculations get boring, let's program!









We start with a specialized program for \( u'=u \), \( u(0)=U_0 \)

Algorithm:

Given \( \Delta t \) (dt) and \( n \)









We start with a specialized program for \( u'=u \), \( u(0)=U_0 \)

Program:

import numpy as np
import sys

dt = float(sys.argv[1])
U0 = 1
T = 4
n = int(T/dt)

t = np.zeros(n+1)
u = np.zeros(n+1)

t[0] = 0
u[0] = U0
for k in range(n):
    t[k+1] = t[k] + dt
    u[k+1] = (1 + dt)*u[k]

# plot u against t









The solution if we plot \( u \) against \( t \)

\( \Delta t = 0.4 \) and \( \Delta t=0.2 \):













The algorithm for the general ODE \( u'=f(u,t) \)

Algorithm:

Given \( \Delta t \) (dt) and \( n \)









Implementation of the general algorithm for \( u'=f(u,t) \)

General function:

def ForwardEuler(f, U0, T, n):
    """Solve u'=f(u,t), u(0)=U0, with n steps until t=T."""
    import numpy as np
    t = np.zeros(n+1)
    u = np.zeros(n+1)  # u[k] is the solution at time t[k]

    u[0] = U0
    t[0] = 0
    dt = T/float(n)

    for k in range(n):
        t[k+1] = t[k] + dt
        u[k+1] = u[k] + dt*f(u[k], t[k])

    return u, t

Magic:

This simple function can solve any ODE (!)









Example on using the function

Mathematical problem:

Solve \( u'=u \), \( u(0)=1 \), for \( t\in [0,4] \), with \( \Delta t = 0.4 \)
Exact solution: \( u(t)=e^t \).

Basic code:

def f(u, t):
    return u

U0 = 1
T = 3
n = 30
u, t = ForwardEuler(f, U0, T, n)

Compare exact and numerical solution:

from scitools.std import plot, exp
u_exact = exp(t)
plot(t, u, 'r-', t, u_exact, 'b-',
     xlabel='t', ylabel='u', legend=('numerical', 'exact'),
     title="Solution of the ODE u'=u, u(0)=1")









Now you can solve any ODE!

Recipe:

Warning:

The Forward Euler method may give very inaccurate solutions if \( \Delta t \) is not sufficiently small. For some problems (like \( u''+u=0 \)) other methods should be used.









Let us make a class instead of a function for solving ODEs

Usage of the class:

method = ForwardEuler(f, dt)
method.set_initial_condition(U0, t0)
u, t = method.solve(T)
plot(t, u)

How?









The code for a class for solving ODEs (part 1)

import numpy as np

class ForwardEuler_v1:
    def __init__(self, f, dt):
        self.f, self.dt = f, dt

    def set_initial_condition(self, U0):
        self.U0 = float(U0)









The code for a class for solving ODEs (part 2)

class ForwardEuler_v1:
    ...
    def solve(self, T):
        """Compute solution for 0 <= t <= T."""
        n = int(round(T/self.dt))  # no of intervals
        self.u = np.zeros(n+1)
        self.t = np.zeros(n+1)
        self.u[0] = float(self.U0)
        self.t[0] = float(0)

        for k in range(self.n):
            self.k = k
            self.t[k+1] = self.t[k] + self.dt
            self.u[k+1] = self.advance()
        return self.u, self.t

    def advance(self):
        """Advance the solution one time step."""
        # Create local variables to get rid of "self." in
        # the numerical formula
        u, dt, f, k, t = self.u, self.dt, self.f, self.k, self.t

        unew = u[k] + dt*f(u[k], t[k])
        return unew









Alternative class code for solving ODEs (part 1)

# Idea: drop dt in the constructor.
# Let the user provide all time points (in solve).

class ForwardEuler:
    def __init__(self, f):

        # test that f is a function
        if not callable(f):
            raise TypeError('f is %s, not a function' % type(f))
        self.f = f

    def set_initial_condition(self, U0):
        self.U0 = float(U0)

    def solve(self, time_points):
        ...









Alternative class code for solving ODEs (part 2)

class ForwardEuler:
    ...
    def solve(self, time_points):
        """Compute u for t values in time_points list."""
        self.t = np.asarray(time_points)
        self.u = np.zeros(len(time_points))

        self.u[0] = self.U0

        for k in range(len(self.t)-1):
            self.k = k
            self.u[k+1] = self.advance()
        return self.u, self.t

    def advance(self):
        """Advance the solution one time step."""
        u, f, k, t = self.u, self.f, self.k, self.t

        dt = t[k+1] - t[k]
        unew = u[k] + dt*f(u[k], t[k])
        return unew









Verifying the class implementation; mathematics

Mathematical problem:

Important result: the numerical method (and most others) will exactly reproduce \( u \) if it is linear in \( t \) (!): $$ \begin{align*} u(t) &= at+b = 0.2t + 3\\ h(t) &= u(t)\\ u'(t) &=0.2 + (u-h(t))^4,\quad u(0)=3,\quad t\in [0,3] \end{align*} $$ This \( u \) should be reproduced to machine precision for any \( \Delta t \).









Verifying the class implementation; implementation

Code:

def test_ForwardEuler_against_linear_solution():
    def f(u, t):
        return 0.2 + (u - h(t))**4

    def h(t):
        return 0.2*t + 3

    solver = ForwardEuler(f)
    solver.set_initial_condition(U0=3)
    dt = 0.4; T = 3; n = int(round(float(T)/dt))
    time_points = np.linspace(0, T, n+1)
    u, t = solver.solve(time_points)
    u_exact = h(t)
    diff = np.abs(u_exact - u).max()
    tol = 1E-14
    success = diff < tol
    assert success









Using a class to hold the right-hand side \( f(u,t) \)

Mathematical problem:

$$ u^{\prime}(t)=\alpha u(t)\left( 1-\frac{u(t)}{R}\right),\quad u(0)=U_0,\quad t\in [0,40]$$

Class for right-hand side \( f(u,t) \):

class Logistic:
    def __init__(self, alpha, R, U0):
        self.alpha, self.R, self.U0 = alpha, float(R), U0

    def __call__(self, u, t):   # f(u,t)
        return self.alpha*u*(1 - u/self.R)

Main program:

problem = Logistic(0.2, 1, 0.1)
time_points = np.linspace(0, 40, 401)
method = ForwardEuler(problem)
method.set_initial_condition(problem.U0)
u, t = method.solve(time_points)









Figure of the solution













Numerical methods for ordinary differential equations

Forward Euler method:

$$ u_{k+1} = u_k + \Delta t\, f(u_k, t_k) $$

4th-order Runge-Kutta method:

$$ u_{k+1} = u_k + {1\over 6}\left( K_1 + 2K_2 + 2K_3 + K_4\right)$$ $$ \begin{align*} K_1 &= \Delta t\,f(u_k, t_k)\\ K_2 &= \Delta t\,f(u_k + {1\over2}K_1, t_k + {1\over2}\Delta t)\\ K_3 &= \Delta t\,f(u_k + {1\over2}K_2, t_k + {1\over2}\Delta t)\\ K_4 &= \Delta t\,f(u_k + K3, t_k + \Delta t) \end{align*} $$

And lots of other methods! How to program a wide collection of methods? Use object-oriented programming!









A superclass for ODE methods

Common tasks for ODE solvers:

Principles:









The superclass code

class ODESolver:
    def __init__(self, f):
        self.f = f

    def advance(self):
        """Advance solution one time step."""
        raise NotImplementedError  # implement in subclass

    def set_initial_condition(self, U0):
        self.U0 = float(U0)

    def solve(self, time_points):
        self.t = np.asarray(time_points)
        self.u = np.zeros(len(self.t))
        # Assume that self.t[0] corresponds to self.U0
        self.u[0] = self.U0

        # Time loop
        for k in range(n-1):
            self.k = k
            self.u[k+1] = self.advance()
        return self.u, self.t

    def advance(self):
        raise NotImplemtedError # to be impl. in subclasses









Implementation of the Forward Euler method

Subclass code:

class ForwardEuler(ODESolver):
    def advance(self):
        u, f, k, t = self.u, self.f, self.k, self.t

        dt = t[k+1] - t[k]
        unew = u[k] + dt*f(u[k], t)
        return unew

Application code for \( u'-u=0 \), \( u(0)=1 \), \( t\in [0,3] \), \( \Delta t=0.1 \):

from ODESolver import ForwardEuler
def test1(u, t):
    return u

method = ForwardEuler(test1)
method.set_initial_condition(U0=1)
u, t = method.solve(time_points=np.linspace(0, 3, 31))
plot(t, u)









The implementation of a Runge-Kutta method

Subclass code:

class RungeKutta4(ODESolver):
    def advance(self):
        u, f, k, t = self.u, self.f, self.k, self.t

        dt = t[k+1] - t[k]
        dt2 = dt/2.0
        K1 = dt*f(u[k], t)
        K2 = dt*f(u[k] + 0.5*K1, t + dt2)
        K3 = dt*f(u[k] + 0.5*K2, t + dt2)
        K4 = dt*f(u[k] + K3, t + dt)
        unew = u[k] + (1/6.0)*(K1 + 2*K2 + 2*K3 + K4)
        return unew

Application code (same as for ForwardEuler):

from ODESolver import RungeKutta4
def test1(u, t):
    return u

method = RungeKutta4(test1)
method.set_initial_condition(U0=1)
u, t = method.solve(time_points=np.linspace(0, 3, 31))
plot(t, u)









The user should be able to check intermediate solutions and terminate the time stepping

def terminate(u, t, step_no):
    eps = 1.0E-6                     # small number
    return abs(u[step_no,0]) < eps   # close enough to zero?









Systems of differential equations (vector ODE)

$$ \begin{align*} u' &= v\\ v' &= -u\\ u(0) &= 1\\ v(0) &= 0 \end{align*} $$













Example on a system of ODEs (vector ODE)

Two ODEs with two unknowns \( u(t) \) and \( v(t) \): $$ \begin{align*} u'(t) &= v(t)\\ v'(t) &= -u(t) \end{align*} $$

Each unknown must have an initial condition, say $$ u(0)=0,\quad v(0)=1 $$

In this case, one can derive the exact solution to be $$ u(t)=\sin (t),\quad v(t)=\cos (t)$$

Systems of ODEs appear frequently in physics, biology, finance, ...









The ODE system that is the final project in the course

Model for spreading of a disease in a population: $$ \begin{align*} S' &= - \beta SI\\ I' &= \beta SI -\nu R\\ R' &= \nu I \end{align*} $$

Initial conditions: $$ \begin{align*} S(0) &= S_0\\ I(0) &= I_0\\ R(0) &= 0 \end{align*} $$









Another example on a system of ODEs (vector ODE)

Second-order ordinary differential equation, for a spring-mass system (from Newton's second law): $$ mu'' + \beta u' + ku = 0,\quad u(0)=U_0,\ u'(0)=0$$

We can rewrite this as a system of two first-order equations, by introducing two new unknowns $$ u^{(0)}(t) \equiv u(t),\quad u^{(1)}(t) \equiv u'(t)$$

The first-order system is then $$ \begin{align*} {d\over dt} u^{(0)}(t) &= u^{(1)}(t)\\ {d\over dt} u^{(1)}(t) &= -m^{-1}\beta u^{(1)} - m^{-1}ku^{(0)} \end{align*} $$

Initial conditions: \( u^{(0)}(0) = U_0 \), \( u^{(1)}(0)=0 \)









Making a flexible toolbox for solving ODEs









Vector notation for systems of ODEs: unknowns and equations

General software for any vector/scalar ODE demands a general mathematical notation. We introduce \( n \) unknowns $$ u^{(0)}(t), u^{(1)}(t), \ldots, u^{(n-1)}(t) $$ in a system of \( n \) ODEs: $$ \begin{align*} {d\over dt}u^{(0)} &= f^{(0)}(u^{(0)}, u^{(1)}, \ldots, u^{(n-1)}, t)\\ {d\over dt}u^{(1)} &= f^{(1)}(u^{(0)}, u^{(1)}, \ldots, u^{(n-1)}, t)\\ \vdots &=& \vdots\\ {d\over dt}u^{(n-1)} &= f^{(n-1)}(u^{(0)}, u^{(1)}, \ldots, u^{(n-1)}, t) \end{align*} $$









Vector notation for systems of ODEs: vectors

We can collect the \( u^{(i)}(t) \) functions and right-hand side functions \( f^{(i)} \) in vectors: $$ u = (u^{(0)}, u^{(1)}, \ldots, u^{(n-1)})$$ $$ f = (f^{(0)}, f^{(1)}, \ldots, f^{(n-1)})$$

The first-order system can then be written $$ u' = f(u, t),\quad u(0) = U_0$$ where \( u \) and \( f \) are vectors and \( U_0 \) is a vector of initial conditions

The magic of this notation:

Observe that the notation makes a scalar ODE and a system look the same, and we can easily make Python code that can handle both cases within the same lines of code (!)









How to make class ODESolver work for systems of ODEs

Forward Euler applied to a system: $$ \underbrace{u_{k+1}}_{\mbox{vector}} = \underbrace{u_k}_{\mbox{vector}} + \Delta t\, \underbrace{f(u_k, t_k)}_{\mbox{vector}} $$

In Python code:

unew = u[k] + dt*f(u[k], t)

where









The adjusted superclass code (part 1)

To make ODESolver work for systems:

class ODESolver:
    def __init__(self, f):
        # Wrap user's f in a new function that always
        # converts list/tuple to array (or let array be array)
        self.f = lambda u, t: np.asarray(f(u, t), float)

    def set_initial_condition(self, U0):
        if isinstance(U0, (float,int)):  # scalar ODE
            self.neq = 1                 # no of equations
            U0 = float(U0)
        else:                            # system of ODEs
            U0 = np.asarray(U0)
            self.neq = U0.size           # no of equations
        self.U0 = U0









The superclass code (part 2)

class ODESolver:
    ...
    def solve(self, time_points, terminate=None):
        if terminate is None:
            terminate = lambda u, t, step_no: False

        self.t = np.asarray(time_points)
        n = self.t.size
        if self.neq == 1:  # scalar ODEs
            self.u = np.zeros(n)
        else:              # systems of ODEs
            self.u = np.zeros((n,self.neq))

        # Assume that self.t[0] corresponds to self.U0
        self.u[0] = self.U0

        # Time loop
        for k in range(n-1):
            self.k = k
            self.u[k+1] = self.advance()
            if terminate(self.u, self.t, self.k+1):
                break  # terminate loop over k
        return self.u[:k+2], self.t[:k+2]

All subclasses from the scalar ODE works for systems as well









Example on how to use the general class hierarchy

Spring-mass system formulated as a system of ODEs:

$$ mu'' + \beta u' + ku = 0,\quad u(0),\ u'(0)\hbox{ known}$$ $$ \begin{align*} u^{(0)} &= u,\quad u^{(1)} = u'\\ u(t) &= (u^{(0)}(t), u^{(1)}(t))\\ f(u,t)&=(u^{(1)}(t), -m^{-1}\beta u^{(1)} - m^{-1}ku^{(0)})\\ u'(t) &= f(u,t) \end{align*} $$

Code defining the right-hand side:

def myf(u, t):
    # u is array with two components u[0] and u[1]:
    return [u[1],
            -beta*u[1]/m - k*u[0]/m]









Alternative implementation of the \( f \) function via a class

Better (no global variables):

class MyF:
    def __init__(self, m, k, beta):
        self.m, self.k, self.beta = m, k, beta

    def __call__(self, u, t):
        m, k, beta = self.m, self.k, self.beta
        return [u[1], -beta*u[1]/m - k*u[0]/m]

Main program:

from ODESolver import ForwardEuler
# initial condition:
U0 = [1.0, 0]
f = MyF(1.0, 1.0, 0.0)   # u'' + u = 0 => u(t)=cos(t)
solver = ForwardEuler(f)
solver.set_initial_condition(U0)

T = 4*pi; dt = pi/20; n = int(round(T/dt))
time_points = np.linspace(0, T, n+1)
u, t = solver.solve(time_points)

# u is an array of [u0,u1] arrays, plot all u0 values:
u0_values = u[:,0]
u0_exact = cos(t)
plot(t, u0_values, 'r-',  t, u0_exact, 'b-')









Throwing a ball; ODE model

Newton's 2nd law for a ball's trajectory through air leads to $$ \begin{align*} {dx\over dt} &= v_x\\ {dv_x\over dt} &= 0\\ {dy\over dt} &= v_y\\ {dv_y\over dt} &= -g \end{align*} $$ Air resistance is neglected but can easily be added!









Throwing a ball; code

Define the right-hand side:

def f(u, t):
    x, vx, y, vy = u
    g = 9.81
    return [vx, 0, vy, -g]

Main program:

from ODESolver import ForwardEuler
# t=0: prescribe x, y, vx, vy
x = y = 0                     # start at the origin
v0 = 5;  theta = 80*pi/180    # velocity magnitude and angle
vx = v0*cos(theta)
vy = v0*sin(theta)
# Initial condition:
U0 = [x, vx, y, vy]

solver= ForwardEuler(f)
solver.set_initial_condition(u0)
time_points = np.linspace(0, 1.2, 101)
u, t = solver.solve(time_points)

# u is an array of [x,vx,y,vy] arrays, plot y vs x:
x = u[:,0];  y = u[:,2]
plot(x, y)









Throwing a ball; results

Comparison of exact and Forward Euler solutions













Logistic growth model; ODE and code overview

Model:

$$ u' = \alpha u(1 - u/{\color{red}R(t)}),\quad u(0)=U_0 $$ \( R \) is the maximum population size, which can vary with changes in the environment over time

Implementation features:









Logistic growth model; class Problem (\( f \))

class Problem:
    def __init__(self, alpha, R, U0, T):
        self.alpha, self.R, self.U0, self.T = alpha, R, U0, T

    def __call__(self, u, t):
        """Return f(u, t)."""
        return self.alpha*u*(1 - u/self.R(t))

    def terminate(self, u, t, step_no):
        """Terminate when u is close to R."""
        tol = self.R*0.01
        return abs(u[step_no] - self.R) < tol

problem = Problem(alpha=0.1, R=500, U0=2, T=130)









Logistic growth model; class Solver

class Solver:
    def __init__(self, problem, dt,
                 method=ODESolver.ForwardEuler):
        self.problem, self.dt = problem, dt
        self.method = method

    def solve(self):
        solver = self.method(self.problem)
        solver.set_initial_condition(self.problem.U0)
        n = int(round(self.problem.T/self.dt))
        t_points = np.linspace(0, self.problem.T, n+1)
        self.u, self.t = solver.solve(t_points,
                                      self.problem.terminate)

    def plot(self):
        plot(self.t, self.u)

problem = Problem(alpha=0.1, U0=2, T=130,
                  R=lambda t: 500 if t < 60 else 100)
solver = Solver(problem, dt=1.)
solver.solve()
solver.plot()
print 'max u:', solver.u.max()









Logistic growth model; results