• Scalar ODE: a single ODE, one unknown function
• Vector ODE or systems of ODEs: several ODEs, several unknown functions

 
$$\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*}$$

• Our methods and software should be applicable to any ODE
• Therefore we need an abstract notation for an arbitrary ODE

 
$$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

• Numerical differentiation: $f'(x)$
• Numerical integration: $\int_a^b f(x)dx$
• Numerical solution of algebraic equations: $f(x)=0$

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$

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})$

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.

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!

 
$$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$

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

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$?

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!

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

 
$$u'=u$$

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

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...

• Exact solution: $u(t)=2e^t$, $u(t_k)=2e^{k\Delta t} = 2(e^{\Delta t})^k$
• Numerical solution: $u_k = 2(1+\Delta t)^k$

When going from $t_k$ to $t_{k+1}$, the solution is amplified by a factor:

• Exact: $u(t_{k+1}) = e^{\Delta t}u(t_k)$
• Numerical: $u_{k+1} = (1+\Delta t)u_k$

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$.

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...

When hand calculations get boring, let's program!

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

• Create arrays t and u of length $n+1$
• Set initial condition: u[0] = $U_0$, t[0]=0
• For $k=0,1,2,\ldots,n-1$:

• t[k+1] = t[k] + dt
• u[k+1] = (1 + dt)*u[k]

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

• Create arrays t and u of length $n+1$
• Create array u to hold $u_k$ and
• Set initial condition: u[0] = $U_0$, t[0]=0
• For $k=0,1,2,\ldots,n-1$:

• u[k+1] = u[k] + dt*f(u[k], t[k]) (the only change!)
• t[k+1] = t[k] + dt

This simple function can solve any ODE (!)

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

• Identify $f(u,t)$ in your ODE
• Make sure you have an initial condition $U_0$
• Implement the $f(u,t)$ formula in a Python function f(u, t)
• Choose $\Delta t$ or no of steps $n$
• Call u, t = ForwardEuler(f, U0, T, n)
• plot(t, u)

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.

• Store $f$, $\Delta t$, and the sequences $u_k$, $t_k$ as attributes
• Split the steps in the ForwardEuler function into four methods:

• the constructor (__init__)
• set_initial_condition for $u(0)=U_0$
• solve for running the numerical time stepping
• advance for isolating the numerical updating formula
  (new numerical methods just need a different advance method, the rest is the same)

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$.

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

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

 
$$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!

• Sometimes a property of the solution determines when to stop the solution process: e.g., when $u < 10^{-7}\approx 0$.
• Extension: solve(time_points, terminate)
• terminate(u, t, step_no) is called at every time step, is user-defined, and returns True when the time stepping should be terminated
• Last computed solution is u[step_no] at time t[step_no]

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

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, ...

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*}$$

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$

• For scalar ODEs we could make one general class hierarchy to solve "all" problems with a range of methods
• Can we easily extend class hierarchy to systems of ODEs?
• Yes!
• The example here can easily be extended to professional code (Odespy)

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*}$$

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

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 (!)

• Recall: ODESolver was written for a scalar ODE
• Now we want it to work for a system $u'=f$, $u(0)=U_0$, where $u$, $f$ and $U_0$ are vectors (arrays)
• What are the problems?

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

• u is a two-dim. array (u[k] is a row)
• f is a function returning an array (all the right-hand sides $f^{(0)},\ldots,f^{(n-1)}$)

• Ensure that f(u,t) returns an array.
  This can be done be a general adjustment in the superclass!
• Inspect $U_0$ to see if it is a number or list/tuple and make corresponding u 1-dim or 2-dim array

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

 
$$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*}$$

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!

• 4 ODEs with 4 unknowns:

• the ball's position $x(t)$, $y(t)$
• the velocity $v_x(t)$, $v_y(t)$

Comparison of exact and Forward Euler solutions

 
$$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

• Class Problem holds "all physics": $\alpha$, $R(t)$, $U_0$, $T$ (end time), $f(u,t)$ in ODE
• Class Solver holds "all numerics": $\Delta t$, solution method; solves the problem and plots the solution
• Solve for $t\in [0,T]$ but terminate when $|u-R| < \hbox{tol}$