<a href="https://colab.research.google.com/github/cohmathonc/biosci670/blob/master/IntroductionComputationalMethods/05_IntroCompMethods_SolvingODEs.ipynb" target="_parent"></a>


```python
import numpy as np
import matplotlib.pylab as plt
import pandas as pd
```

# Numerical Solution of ODEs

Consider a first order differential equation of the form

$$y'=\frac{dy}{dx}=f(x, y)\, , \qquad y(x=0)=y_0 \tag{1}$$

We want to find the function $y(x)$ that describes the evolution of this quantity in function of $x$.
Let's assume the value of $y(x_0)$ at a specific point $x=x_0$ is known;  $y(x=x_0)$ is called the *initial value*.

The task of finding function $y(x)$, given its derivative and initial value as in Eq. (1), is called an *initial value problem*.

Recall the definition of the continuous derivative of an analytic function $f(x)$:

$$y'(x_0)=\lim_{\Delta x\rightarrow 0}\frac{y(x_0+\Delta x)-y(x_0)}{\Delta x} =f(x_0, y)\; . \tag{2}$$


## Numerical Integration Methods

We approximate Eq (2) by finite differences:

$$y'(x)\approx \frac{y(x_0+\Delta x)-y(x_0)}{\Delta x} =f(x_0, y)\; , \tag{3}$$

which allows us to express $y$ at a neighbouring point $x=x_0+\Delta x$ in terms of its value at the current position $x_0$ and its derivative at this point:

$$y(x_0+\Delta x) = y(x_0) + y'(x)\, \Delta x = y(x_0) + f(x_0, y(x_0))\, \Delta x  \tag{4} \, .$$

We can make the following geometric analogy:

Consider the problem of calculating the shape of a curve $y(x)$ which starts at a given point $y(x_0)$and satisfies the differential equation (1). 
The differential equation can be thought of as a formula that describes the slope of the tangent line to the curve $y$ at any point $x_i$, provided that the position of that point on the curve $y(x_i)$ is known.
Given a starting point $x_0$, we can compute the next point on the curve $y(x_0 + \Delta x)=y(x_1)$ using Eq. (4). This procedure can be repeated until the end of the domain is reached.


In this course we are mainly interested in the **time**-evolution of biological systems.
Therefore, for our purpose, the variable $x$ in (1) and (2) will be time $t$ so that (4) becomes:

$$y(t_0+\Delta t) = y(t_0) + y'(t)\, \Delta t = y(t_0) + f(t_0, y(t_0))\, \Delta t  \tag{5} \, .$$



### Euler Method (explicit)

Assume that we wish to approximate the unknown function $y(t)$ in the interval $[t_{\text{min}}, t_\text{max}]$ which is discretized into $N$ subintervals with $t_{\text{min}}=t_0 < t_1 < t_2 < \ldots <t_{N-1} < t_{N} =t_{\text{max}}$. 
We define the *step size* between subsequent points as $\Delta_i = t_{i+1}-t_i$,and
write Eq. (4) as a *recursive* scheme:

$$ y_{i+1} = y_i + \Delta_i\, f(t_i, y_i) \tag{5}$$

This is called the **Euler method**. 
The Euler method is an *explicit* integration method because any $y_{i+1}$ is an explicit function of the 'previous' $y_i$ for $i\leq N-1$:

\begin{align}
y_1 &= y_0 + \Delta_0 \, f(t_0, y_0) \\
y_2 &= y_1 + \Delta_1 \, f(t_1, y_1) \\
& \dots \\
y_{N-1} &= y_{N-2} + \Delta_{N-2} \, f(t_{N-2}, y_{N-2}) \\
y_{N} &= y_{N-1} + \Delta_{N-1} \, f(t_{N-1}, y_{N-1}) \\
\end{align}

If the evaluation points are spaced equally, $\Delta_i = \Delta = \left| t_{\text{max}}-t_\text{min}\right|/N$.


The following example illustrates how this method approximates $y(t)$.


```python
#=================================================
# Approximate function using Euler's method
#
# Explore effect of varying 'delta_t'
#=================================================
delta_t = 0.2
t_0 = 0
y_0 = 1


# function defs & preparation
from cycler import cycler
colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728',
          '#9467bd', '#8c564b', '#e377c2', '#7f7f7f',
          '#bcbd22', '#17becf']
color_cyler = cycler(color=colors)

def f_y(y_0, t):
    return y_0*np.exp(2*t)
    
def f_yp(y, t):
    return 2*y

def plot_deriv(ax, x, y, yp, delta_x=0.1, **kwargs):
    x_plot = np.array([x, x+delta_x])
    y_plot = np.array([y, y+yp*delta_x])
    ax.plot(x_plot[1], y_plot[1], '.k')
    ax.plot(x_plot, y_plot, linestyle = '--', **kwargs)
    
cm = plt.get_cmap('viridis_r')
n_steps = int(1/delta_t)
cgen = (cm(1.*i/n_steps) for i in range(n_steps))
color_cyler = cycler(color=cgen)

# plot (unknown) function
t = np.arange(t_0,1,0.05)
y = f_y(y_0, t)
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111) 
ax.plot(t, y, ':', label='y (unknown)')

# plot approximation
t_i = t_0
y_i = y_0
ax.plot(t_i, y_i, 'sk', label="y(t=%.2f)"%t_i)
for i, kwargs in zip(range(1,n_steps), color_cyler):
    plot_deriv(ax, t_i, y_i, f_yp(y_i, t_i), delta_x=delta_t, **kwargs)
    t_next = t_i + delta_t
    y_next = y_i + delta_t*f_yp(y_i, t_i)
    ax.plot(t_next, y_next, marker='s', **kwargs)
    t_i = t_next
    y_i = y_next
ax.legend()   
```



---
**Exercise (1):**

1. Write a function that computes: 
  $$\frac{d y(t)}{d t} = f(t,y)= y\, t\, .$$

2. Write a function that implements the Euler Method.
  This function should accept the following arguments:

  - *f* : a function that computes $\frac{dy}{dt}$ in function of $y$ and $t$, i.e. $f(t,y)$ above. 
  - *t*: array of points (evaluation time points) at which y(t) should be estimated
  - *y_0*: initial value at $t=0$ 
  
---



```python
def dydt(t, y):
  """  
  Args:
    - t: scalar, usually time
    - y: scalar, function value
    
  Returns:
    - scalar, derivative dy/dt of (unknown) function y(t)
  """
  return y*t

def solve_euler(f, t, y_0):
  """
  Uses explicit euler method to solve ODE: y'=f(t, y) 
  with initial value y(t_0)=y_0.
  
  Args:
    - f: Function object that computes y', expected signature f(t, y), 
         where t is evaluation point and y is function value
    - t: array of evaluation points
    - y_0: initial value, i.e. y(t[0])
  
  Returns:
    - array containing approximated function values y(t[i]) 
  """
  y = np.zeros_like(t)
  y[0] = y_0
  for i in range(0, len(t)-1):
    y[i+1] = y[i] + (t[i+1]-t[i]) * f(t[i],y[i])
  return y

```


---
**Exercise (2):**

1. Use this function to solve 
$$\frac{dy}{dt}=r\, y \qquad y(t_0=0)=y_0=0$$
in the interval $t\in [0, 10]$, with constant $r=0.5$.

2. Compare your approximate solutions to the analytic solution for different choices of integration step size:
$$y(t) = y_0\,  e^{r\, t}$$

  What step size do you need for a good approximation?
  
3. Solve above equation for $r=-2.5$ and compare to analytic solution. What  do you observe at large step sizes $\Delta t \gtrapprox 0.5$.



---

*Note regarding function arguments*:

The ODE in above example describes an *autonomous system*:
$$\frac{d y }{dt} = f(y(t), t) = f(y(t))$$
that does not explicitly depend on the independent variable, time in this case.

The computational function `f()` for computing the value $\frac{dy}{dt}$ for a given $y$ therefore does not require $t$ as an argument.
However, the implementation of the ODE solver function should be sufficiently generic so that it can be applied to autonomous as well as non-autonomous (with time dependency) systems.

Standard ODE solvers therefore expect the function `f()` that computes the derivative  to accept two arguments, the current (estimated) value $y_i$ and the current time $t$. 
Different conventions exist for the order in which those arguments are expected, either `f(y, t)` or `f(t, y)` respectively.
The order and type of the arguments of a function is called the *signature* of the function.
We use the second option `f(t, y)` throughout this notebook, i.e. a function signature that has time as first argument, and function values as second argument.
Also, we will include both arguments even when describing an autonomous system.



```python
# define functions
def dydt_1(t, y):
  return 0.5*y

# analytic:
def fun_an_1(t, y_0):
  return y_0*np.exp(0.5*t)


# define functions
def dydt_2(t, y):
  return -2.5*y

# analytic:
def fun_an_2(t, y_0):
  return y_0*np.exp(-2.5*t)


# domain bounds
t_0 = 0
t_N = 10

# initial value
y_0 = 1

# create time data for numeric approximation
n_steps = 15
t = np.linspace(t_0, t_N, n_steps)
step_size = (t_N-t_0)/n_steps
#print("step size: ",step_size)

# solve ODE
y_1 = solve_euler(f=dydt_1, t=t, y_0=y_0)
y_2 = solve_euler(f=dydt_2, t=t, y_0=y_0)

# create data for analytic function on fine grid
t_an = np.linspace(t_0, t_N, 100)
y_an_1 = fun_an_1(t_an, y_0)
y_an_2 = fun_an_2(t_an, y_0)

# plot
fig, axes= plt.subplots(1, 2, figsize=plt.figaspect(0.3))
axes[0].plot(t, y_1, label='numeric approximation of $y(t)$')
axes[0].plot(t_an, y_an_1, label='y(t)')
axes[0].set_ylabel("y(t)")
axes[0].set_xlabel("t")
axes[0].set_title("function (1): $r=0.5$; stepsize $\Delta t = %.2f$"%step_size)
axes[0].legend()

axes[1].plot(t, y_2, label='numeric approximation of $y(t)$')
axes[1].plot(t_an, y_an_2, label='y(t)')
axes[1].set_ylabel("y(t)")
axes[1].set_xlabel("t")
axes[1].set_title("function (2): $r=-2.5$; stepsize $\Delta t = %.2f$"%step_size)
axes[1].legend()

plt.show()
```

In the first of above examples ($r=0.5$), we noticed that the numeric approximation underestimates the values of the actual function $y(t)$ for stepsizes $\Delta t \gtrapprox 0.01$.
In the second example ($r=-2.5$), we saw that for step sizes $\Delta t \gtrapprox 0.5$, the numeric approximation may become *unstable*, assuming very large values or oscillating around the actual function value.

To ensure convergence with the explicit Euler method, very small stepsizes must be used. Oscillating function values can be numeric artefacts and a sign of numeric instability!




#### Error and Convergence

The Euler method is an iterative method. In each iteration step, it estimates the function value $y_i$ using a finite difference approximation, Eq. (4), that involves previous (estimated) values $y_{i-1}$.
As we have discussed in the notebook on [numeric derivatives](https://github.com/cohmathonc/biosci670/blob/master/IntroductionComputationalMethods/03_IntroCompMethods_NumericalDifferentiation.ipynb), this approximation involves a *truncation error* because it neglects contributions of higher function derivatives. 

The truncation error that is made in iteration step $i+1$, when there was no error in step $i$, is called *local truncation error* (LTE).
To compute the LTE, we compare the Taylor expansion of $y(t+\Delta t)$ to the iterative formula of the Euler method, Eq (5):
\begin{align}
 y(t_i+\Delta t) &= y(t_i) + \Delta t\, y'(t_i) + \frac{1}{2}\Delta t^2\, y''(\tau) + \mathcal{O}(\Delta t^3) \tag{6a}\\
 y_{i+1} &= y_i + \Delta t \underbrace{f(t_i, y_i)}_{y'} \tag{6b}
\end{align}

where we assume that no error has been made in previous steps, i.e. that the estimate $y_i$ is identical to the actual function value $y(t_i)$. The difference between (6a) and (6b) corresponds to the LTE:
$$\text{LTE}:\qquad \epsilon^{\text{local}}_i = \frac{1}{2}\Delta t^2\, y''(\tau) + \mathcal{O}(\Delta t^3) \tag{7}$$

In an iterative method, errors accumulate with each additional step. The cumulative effect of these errors is called the *global truncation error* (GTE). The GTE is the error at a fixed time $t$ and after however many steps the method needs to take to reach this time from an initial time point $t_0$:
$$\text{GTE}:\qquad \epsilon^{\text{global}}_i = y(t_i) - y_i$$

For a given integration interval $[t_0, t_\text{max}]$, the number of steps is proportional to $\frac{1}{\Delta t}$, and it can be shown that the GTE is proportional to $\Delta t$. The Euler method is therefore a first order method.



---
**Exercise (3):**

Investigate the convergence behavior of the Euler method by computing the approximation error
$$\epsilon_i = y(t_i) - y_i$$
- at the very first iteration ($i=1$) step to obtain an estimate of the LTE (first step, no accumulation of errors) and 
- at some later evaluation step to obtain an estimate of the GTE.

Repeat this computation for multiple values of stepsize $\Delta t$, halfing the stepsize each iteration step.
Plot the error terms in function of stepsize $\Delta t$ and compute the convergence order.

---




```python
# define functions for common analysis steps

def compute_error(f_yt, solver, solution, t_0, t_N, y_0):
  """
  Computes cumulative error between numeric approximation and analytic solution
  to ODE.
  
  Args:
    - f_yt: Function object that computes y', expected signature f(t, y), 
            where t is evaluation point and y is function value
    - solver: ODE solver, expected signature: solver(f_yt, t, y_0)
            where t is array of time points.
    - solution: Function object that computes value y(t) for analytic solution.
            Expected signature: solution(t, y_0)
    - t_0, t_N: endpoints of integration interval
    - y_0: initial value, i.e. y(t[0])
    
  Returns:
    - pandas DataFrame containing 'lte' and 'gte' for different step sizes 'delta_t'
  """
  df = pd.DataFrame()
  for i in range(0,20):
      n_steps = 1 * (2**i)
      delta_t = (t_N - t_0)/n_steps
      t = np.linspace(t_0, t_N, n_steps)
      y_approx = solver(f_yt, t, y_0)
      y_exact  = solution(t, y_0)
      error_cum = np.sum(np.abs(y_exact - y_approx))/n_steps
      # keep record of parameters
      df.loc[i, 'n_steps'] = n_steps
      df.loc[i, 'delta_t'] = delta_t
      if n_steps > 2:
          ref_step = n_steps - 1  
          error_gte = np.abs(y_exact[ref_step] - y_approx[ref_step])  
          df.loc[i, 'gte'] = error_gte  
      if n_steps > 1:
          error_lte = np.abs(y_exact[1] - y_approx[1])
          df.loc[i, 'lte'] = error_lte
  return df


def plot_error(df, method):
  """
  Plots cumulative error vs step size in log-log scale.
  
  Args:
    - df: pandas dataframe object with columns 'delta_t', 'lte', 'gte'
    - method: the name of the method that is being plotted
  """
  fig = plt.figure(figsize=plt.figaspect(0.5))
  ax = fig.add_subplot(111) 

  ax.loglog(df.delta_t.values, df.gte.values, marker='o', 
            linestyle='-', label='GTE')
  ax.loglog(df.delta_t.values, df.lte.values, marker='o', 
            linestyle='-', label='LTE')
  ax.set_xlabel('$\Delta_i$')
  ax.set_ylabel("absolute error")
  ax.set_title('Approximation Error in %s'%method)
  ax.legend()
  plt.show()
  
  
def compute_convergence_order(eps):
  """
  Computes ratio of errors steps: ratio(n) = error(n-1)/error(n)
  """
  eps_shifted = np.roll(eps, shift=1)
  eps_shifted[0] = np.nan
  eps_ratio = eps_shifted / eps
  convergence_order = np.log2(eps_ratio)
  return convergence_order
```


```python
# ODE
def dydt(t, y):
  return y

# solution analytic
def fun_an(t, y_0):
  return y_0*np.exp(t)

# domain bounds
t_0 = 0
t_N = 10

# initial value
y_0 = 1

# compute error for different number of integration steps
df_euler = compute_error(dydt, solve_euler, fun_an, t_0, t_N, y_0)

# plot
plot_error(df_euler, 'euler method')

# convergence order
print("Convergence order of explicit Euler method (LTE) ", 
     compute_convergence_order(df_euler.lte.values))
print("Convergence order of explicit Euler method (GTE) ", 
     compute_convergence_order(df_euler.gte.values))
```

### Euler Method (implicit)

Instead of formulating the Euler method as an *explicit* iterative scheme, it can also be formulated as an *implicit* method:
$$ y_{i+1} = y_i + \Delta_i\, f(x_{i+1}, y_{i+1}) \tag{8}$$

 where the new appproximated value $y_{i+1}$ appears on both sides of the equation.
 
 This method is more stable than the explicit Euler method, but computationally more expensive because an algebraic equation needs to be solved in each time step for the unknown $y_{i+1}$.
 Like the explicit Euler method, the implicit method is of first order.
 


```python
def solve_euler_implicit(f, t, y_0):
  """
  Uses implicit euler method to solve ODE: y'=f(t, y) 
  with initial value y(t_0)=y_0.
  
  Args:
    - f: Function object that computes y', expected signature f(t, y), 
         where t is evaluation point and y is function value
    - t: array of evaluation points
    - y_0: initial value, i.e. y(t[0])
  
  Returns:
    - array containing approximated function values y(t[i]) 
  """
  y = np.zeros_like(t)
  y[0] = y_0
  for i in range(0, len(t)-1):
    delta_t = t[i+1] - t[i]
    yp = f(t[i+1], y[i]) / (1+delta_t)
    y[i+1] = y[i] + delta_t * yp
  return y
```


```python
# compute error for different number of integration steps
df_euler_implicit = compute_error(dydt, solve_euler_implicit, 
                                  fun_an, t_0, t_N, y_0)

# plot
plot_error(df_euler_implicit, 'euler method (implicit)')

# convergence order
print("Convergence order of implicit Euler method (LTE) ", 
     compute_convergence_order(df_euler_implicit.lte.values))
# convergence order
print("Convergence order of implicit Euler method (GTE): ", 
     compute_convergence_order(df_euler_implicit.gte.values))
```

### Modified Euler Method

The Euler method tries to predict the next point on the curve $y(t_{i+1})$ from the slope of the curve at an earlier evaluation point. If the curve is convex, the estimated value $y(t_{i+1})$ will be lower than the actual value at $t_{i+1}$, if it is concave, the estimated value wil be higher than the actual value.

Instead of using the slope of $y(t_i)$ at $t_i$ or $t_{i+1}$, the *modified Euler methods* try to use better approximations of the slope.
The so-called *Heun's method* follows an approach similar to the *trapezoidal rule*:
It uses an 'optimal slope' that would connect the points $y(t_i)$ and $y(t_{i+1})$ and which equals the average of the slopes at $y(t_i)$ and $y(t_{i+1})$:
$$\text{average slope} = \frac{\Delta y}{\Delta t}=\frac{y_{i+1}-y_i}{t_{i+1}-t_i}=\frac{1}{2}\left(f(t_i, y_i) + f(t_{i+1}, y_{i+1})\right)$$

And therefore 

$$ y_{i+1} = y_i + \text{average slope} \cdot \Delta t = y_i + \frac{1}{2}\Delta t\left(f(t_i, y_i) + f(t_{i+1}, \tilde{y}_{i+1})\right) \, ,$$

with the intermediate value 

$$\tilde{y}_{i+1}= y_i + \Delta t\, f(t_i, y_i)$$ 




A similar approach, based on the *midpoint rule* exists.
Both methods are of second order accuracy.



```python
def solve_euler_modified(f, t, y_0):
  """
  Uses explicit euler method to solve ODE: y'=f(t, y) 
  with initial value y(t_0)=y_0.
  
  Args:
    - f: Function object that computes y', expected signature f(t, y), 
         where t is evaluation point and y is function value
    - t: array of evaluation points
    - y_0: initial value, i.e. y(t[0])
  
  Returns:
    - array containing approximated function values y(t[i]) 
  """
  
  y = np.zeros_like(t)
  y[0] = y_0
  for i in range(0, len(t)-1):
    y_int  = y[i] + (t[i+1]-t[i]) * f(t[i],y[i])
    y[i+1] = y[i] + 0.5*(t[i+1]-t[i]) * ( f(t[i],y[i]) + f(t[i+1],y_int) )
  return y
```


```python
# compute error for different number of integration steps
df_euler_mod = compute_error(dydt, solve_euler_modified, 
                                  fun_an, t_0, t_N, y_0)

# plot
plot_error(df_euler_mod, 'modified euler method')

# convergence order
print("Convergence order of modified Euler method (GTE): ", 
     compute_convergence_order(df_euler_mod.gte.values))
```

### Runge-Kutta Methods

The Euler and modified Euler methods above are specific cases of a  familiy of iterative methods used for approximating solutions to ODEs.
These methods are called [*Runge-Kutta*](https://en.wikipedia.org/wiki/Runge–Kutta_methods) (RK) methods.
Explicit and implicit Euler methods are *one stage* RK methods. The modified Euler methods are *two stage* RK methods.
In general, RK methods use the weighted average of function values evaluated at multiple increments (*stages*) to approximate the next value $y(t_{i+1})$.
The 'classical' Runge Kutta Method is RK4 which employs four stages for approximation.

\begin{align}
y_{i+1} &= y_i + \frac{1}{6} \left( k_1 + 2 k_2, +2k_3 + k_4\right)\\
t_{i+1} &= t_i + \Delta t
\end{align}

where

\begin{align}
k_1 &= \Delta t \, f(t_i, y_i)\\
k_2 &= \Delta t \, f(t_i+\frac{\Delta t}{2}, y_i + \frac{k_1}{2})\\
k_3 &= \Delta t \, f(t_i+\frac{\Delta t}{2}, y_i + \frac{k_2}{2})\\
k_4 &= \Delta t \, f(t_i+\Delta t, y_i + k_3)\\
\end{align}






```python
def solve_RK4(f, t, y_0):
  """
  Uses explicit euler method to solve ODE: y'=f(t, y) 
  with initial value y(t_0)=y_0.
  
  Args:
    - f: Function object that computes y', expected signature f(t, y), 
         where t is evaluation point and y is function value
    - t: array of evaluation points
    - y_0: initial value, i.e. y(t[0])
  
  Returns:
    - array containing approximated function values y(t[i]) 
  """
  
  y = np.zeros_like(t)
  y[0] = y_0
  for i in range(0, len(t)-1):
    delta_t = t[i+1]-t[i]
    k_1  = delta_t * f(t[i], y[i])
    k_2  = delta_t * f(t[i] + delta_t/2, y[i] + k_1/2)
    k_3  = delta_t * f(t[i] + delta_t/2, y[i] + k_2/2)
    k_4  = delta_t * f(t[i] + delta_t, y[i] + k_3)
    y[i+1] = y[i] + 1/6*(k_1 + 2*k_2 + 2*k_3 + k_4)
  return y
```


```python
# compute error for different number of integration steps
df_RK4 = compute_error(dydt, solve_RK4, 
                                  fun_an, t_0, t_N, y_0)

# plot
plot_error(df_RK4, 'RK4 method')

# convergence order
print("Convergence order of RK4 method (GTE): ", 
     compute_convergence_order(df_RK4.gte.values))
```

## Adaptive Stepsize

In all examples above, we used a constant step size $\Delta t$ for integrating $f(t,y)$.
However, this may be inefficient.

An *adaptive* integration method finds the appropriate $\Delta t$ to meet a given precision target. In Eq. (7) we introduced the local truncation error LTE. Suppose, you evaluate an estimate for $y_{i+1}$ twice, using a method with n-th order accuracy (LTE $\mathcal{O}(\Delta t^{n+1})$) yielding estimate $y_{i+1}^{(n)}$, and another method with $(n+1)$-th order accuracy (LTE $\mathcal{O}(\Delta t^{n+2})$) , yielding estimate $y_{i+1}^{(n+1)}$.
The error $\epsilon$ between $(n+1)$ and $n$ order methods can then be estimated as:
$$\epsilon = \left| y_{i+1}^{(n+1)} - y_{i+1}^{(n)}\right| = \mathcal{O}(\Delta t^{n+1})$$

Suppose we are working with 4th- and 5th- order Runge Kutta Methods, i.e. $n=4$ and would like to find the optimal stepsize $\Delta \tau$ so that  $\epsilon \leq \delta$ for some *desired accuracy* of $\delta$.
From above we see that the error $\epsilon$ scales with $\Delta t^5$.
Therefore
$$\frac{\delta}{\epsilon} = \left(\frac{\Delta\tau}{\Delta t}\right)^{5}$$
and solving for $\Delta\tau$
$$\Delta\tau = \Delta t\,  \left(\frac{\delta}{\epsilon}\right)^{1/5}$$

If $\epsilon < \delta$, we accept the current estimate $y_{i+1}^{(5)}$ and continue to the next time step, possibly increasing the step size for the next step.
If $\epsilon > \delta$, the current estimate is rejected and we recompute the failed time step using a smaller value of $\Delta \tau$.
You can read more details about this approach in the [numerical recipes](https://aip.scitation.org/doi/pdf/10.1063/1.4823060). 

We will not implement this ourselves, but the ODE integration functions provided by `scipy.integrate` use adaptive stepsizes for integration. See examples below.

## Comparison: Runge Kutta 1, 2, 4 stage

Here we compare the different integration methods discussed previously by applying them to the ODEs of exercise (2):
- 1 stage Runge Kutta: explicit & implicit Euler
- 2 stage Runge Kutta: explicit modified Euler
- 4 stage Runge Kutta: 'RK4'

We also include two solvers of the `scipy.integrate` package in this comparison:
- [`scipy.integrate.odeint`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html#scipy.integrate.odeint), based on [ODEPACK LSODA](https://computation.llnl.gov/casc/odepack/) which automatically chooses between methods for *stiff* and *non-stiff* systems.
- [`scipy.integrate.RK45`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.RK45.html#scipy.integrate.RK45), an explicit *adaptive* Runge-Kutta method.

We call latter using the standard interface [`scipy.integrate.solve_ivp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp) which provides access to multiple ode-integration methods.
This interface returns not only the estimated values of $y_i$ but a results *object* that contains various types of information. You can inspect the contents of this object with `dir()`.

Note that the interfaces `odeint` and `solve_ivp` expect the arguments of the ODE-defining function to be ordered differently.
All examples in this notebook use the function signature $f(t, y)$ as expected by `solve_ivp`.


```python
# we also include two solvers from the scipy package
from scipy.integrate import solve_ivp   
from scipy.integrate import odeint

# Function (1) 
def dydt_1(t, y):
  return 0.5*y

def fun_an_1(t, y_0):
  return y_0*np.exp(0.5*t)


# Function (2)
def dydt_2(t, y):
  return -2.5*y

def fun_an_2(t, y_0):
  return y_0*np.exp(-2.5*t)


# integration bounds
t_0 = 0
t_N = 10

# initial value
y_0 = 1

# create time data for numeric approximation
n_steps = 30
t = np.linspace(t_0, t_N, n_steps)
step_size = (t_N-t_0)/n_steps
#print("step size: ",step_size)

# create data for analytic function on fine grid
t_an = np.linspace(t_0, t_N, 100)
y_an_f1 = fun_an_1(t_an, y_0)
y_an_f2 = fun_an_2(t_an, y_0)

# solve ODEs
## f1
y_euler_f1 = solve_euler(f=dydt_1, t=t, y_0=y_0)
y_euler_imp_f1 = solve_euler_implicit(f=dydt_1, t=t, y_0=y_0)
y_euler_mod_f1 = solve_euler_modified(f=dydt_1, t=t, y_0=y_0)
y_RK4_f1 = solve_RK4(f=dydt_1, t=t, y_0=y_0)
### scipy solve_ivp
RK45_f1 = solve_ivp(dydt_1, [t_0, t_N], [y_0], 
                        method='RK45')
y_RK45_f1 = RK45_f1.y.flatten()
t_RK45_f1 = RK45_f1.t
### scipy odeint
y_odeint_f1 = odeint(dydt_1, y_0, t, tfirst=True)
y_odeint_f1 = np.array(y_odeint_f1).flatten()
## f2
y_euler_f2 = solve_euler(f=dydt_2, t=t, y_0=y_0)
y_euler_imp_f2 = solve_euler_implicit(f=dydt_2, t=t, y_0=y_0)
y_euler_mod_f2 = solve_euler_modified(f=dydt_2, t=t, y_0=y_0)
y_RK4_f2 = solve_RK4(f=dydt_2, t=t, y_0=y_0)
### scipy solve_ivp
RK45_f2 = solve_ivp(dydt_2, [t_0, t_N], [y_0], 
                        method='RK45')
y_RK45_f2 = RK45_f2.y.flatten()
t_RK45_f2 = RK45_f2.t
### scipy odeint
y_odeint_f2 = odeint(dydt_2, y_0, t, tfirst=True)
y_odeint_f2 = np.array(y_odeint_f2).flatten()


# plot
fig, axes= plt.subplots(1, 2, figsize=plt.figaspect(0.3))
axes[0].plot(t, y_euler_f1, label='Euler explicit')
axes[0].plot(t, y_euler_imp_f1, label='Euler implicit')
axes[0].plot(t, y_euler_mod_f1, label='Euler modified')
axes[0].plot(t, y_RK4_f1, label='RK4')
axes[0].plot(t_RK45_f1, y_RK45_f1, label='scipy RK45 -- adaptive')
axes[0].plot(t,y_odeint_f1, label='scipy odeint')
axes[0].plot(t_an, y_an_f1, label='analytic $y(t)=y_0\, e^{0.5\,t}$')
axes[0].set_title("function (1): $r=0.5$; stepsize $\Delta t = %.2f$"%step_size)
axes[0].set_ylabel("y(t)")
axes[0].set_xlabel("t")
axes[0].legend()
axes[1].plot(t, y_euler_f2, label='Euler explicit')
axes[1].plot(t, y_euler_imp_f2, label='Euler implicit')
axes[1].plot(t, y_euler_mod_f2, label='Euler modified')
axes[1].plot(t, y_RK4_f2, label='RK4')
axes[1].plot(t_RK45_f2,y_RK45_f2, label='scipy RK45 -- adaptive')
axes[1].plot(t,y_odeint_f2, label='scipy odeint')
axes[1].plot(t_an, y_an_f2, label='analytic $y(t)=y_0\, e^{-2.3\,t}$')
axes[1].set_title("function (2): $r=-2.5$; stepsize $\Delta t = %.2f$"%step_size)
axes[1].set_ylabel("y(t)")
axes[1].set_xlabel("t")
axes[1].legend()
plt.show()
```

## Solving Systems of ODEs

The same methods discussed before can also be used to solve systems of ODE's.

Consider the following example of a system of two ODEs:
\begin{align}
\frac{dx}{dt}&=-x+x^3 &=f(t, x) \\
\frac{dy}{dt}&=-2y & = g(t, y)
\end{align}

For numerical implementation, we first need to write a function that returns the values $f(t, x)$ and $g(t, y)$ for given $x, y, t$.
Ideally, we would like this function to have the same signature (i.e. type and order of arguments) as our previous functions, i.e. one argument for function values, and one for time.
This can be achieved by 'vectorizing' these functions:

Instead of accepting only a single scalar as input for the  current function value `y`, we permit an n-dimensional array as input argument. 
Likewise, instead of returning only the derivative $y'$, we return an array of derivative values, one for each ODE in our system of ODEs.


---
**Exercise (4):**

1. 'Vectorize' the function `dydt` from the beginning of this notebook and implement above system of ODEs.

2. Modify the implementation of the explicit Euler method from before to work with this vectorized function. 

3. Test your implementation using different combinations of initial values $x_0$, $y_0$.

---






```python
def dydt_ndim(t, y):
  """  
  Args:
    - y: array of function values
    - t: scalar, time
    
  Returns:
    - array of derivatives dy/dt of (unknown) function y(t)
  """
  #print("in: ", y)
  y_out = np.zeros_like(y)
  y_out[0] = -y[0] + np.power(y[0],3)
  y_out[1] = -2*y[1]
  #print("out: ", y_out)
  return y_out

def solve_euler_ndim(f, t, y_0):
  """
  Uses explicit euler method to solve ODE: y'=f(y, t) 
  with initial value y(t_0)=y_0.
  
  Args:
    - f: Function object that computes y', expected signature f(y, t), 
         where t is evaluation point and y is function value
    - t: array of evaluation points
    - y_0: array of initial values, i.e. y(t[0])
  
  Returns:
    - array containing approximated function values y(t[i]) 
  """
  y = np.zeros((len(y_0),len(t)))
  y[:,0] = y_0
  for i in range(0, len(t)-1):
    y[:,i+1] = y[:,i] + (t[i+1]-t[i]) * f(t[i], y[:,i])
    #print(y[:,i+1])
  return y


```


```python
# domain bounds
t_0 = 0
t_N = 10

# initial value
y_0 = [0.5,1]

# create time data for numeric approximation
n_steps = 1000
t = np.linspace(t_0, t_N, n_steps)
step_size = (t_N-t_0)/n_steps
print("step size: ",step_size)

# solve ODE
y_euler_ndim = solve_euler_ndim(f=dydt_ndim, t=t, y_0=y_0)

# plot
fig = plt.figure(figsize=plt.figaspect(0.5))
ax = fig.add_subplot(111) 
ax.plot(t, y_euler_ndim[0,:], label='x(t)')
ax.plot(t, y_euler_ndim[1,:], label='y(t)')
ax.set_xlabel("t")
plt.legend()
plt.show()
```

# Useful Python / Scipy Functions


```python
from scipy.integrate import odeint

# define functions
def dydt(y, t): # note that function signature differs from our convention 
  return -2.3*y

# domain bounds
t_0 = 0
t_N = 10

# initial value
y_0 = 1

# create time data for numeric approximation
n_steps = 10
t = np.linspace(t_0, t_N, n_steps)

# expects dy/dx to have calling signature: fun(y, t) !
sol = odeint(dydt, y_0, t)
print("output: ", sol)
sol = np.array(sol).flatten()      # reshape output format from odeint
print("'flattened' output: ", sol)
```


```python
from scipy.integrate import solve_ivp

# define functions
def dydt(t, y):   
  return -2.3*y

# domain bounds
t_0 = 0
t_N = 10

# initial value
y_0 = 1

# automatic time-step size: initial, min, max step size can be defined
# expects dy/dx to have calling signature: fun(t, y) !
sol = solve_ivp(dydt, [t_0, t_N], [y_0])   # supports different integration methods, default is 'RK45'
print("results object: ",sol)

```

## Exercises

- In [this](https://github.com/cohmathonc/biosci670/blob/master/IntroductionComputationalMethods/exercises/07_ODEs.ipynb) exercise you use the explicit Euler method (and optionally also other methods) to solve an initial value problem.

###### About 
This notebook is part of the *biosci670* course on *Mathematical Modeling and Methods for Biomedical Science*.
See https://github.com/cohmathonc/biosci670 for more information and material.

04a1661acff7c68181cc3b4da146847f5ddc72a8
ipynb
Jupyter Notebook
IntroductionComputationalMethods/05_IntroCompMethods_SolvingODEs.ipynb
sbranciamore/biosci670coh
92e91aa1e0266b54b0831537b857796eaf0b949c
CC-BY-4.0
2019-02-19T21:20:42.000Z
2021-01-12T20:25:42.000Z
IntroductionComputationalMethods/05_IntroCompMethods_SolvingODEs.ipynb
danielabler/biosci670
a9fc3e3062cc7ee4f2da3c26a4c5fc22c5cb8cb4
CC-BY-4.0
IntroductionComputationalMethods/05_IntroCompMethods_SolvingODEs.ipynb
danielabler/biosci670
a9fc3e3062cc7ee4f2da3c26a4c5fc22c5cb8cb4
CC-BY-4.0
2019-05-15T17:33:41.000Z
2021-02-22T18:11:37.000Z
Qwen/Qwen-72B
1. YES
2. YES