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


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

# Numeric Solution of Logistic Growth Problem

$$\frac{d\, N}{d\, t}=r\, N \left(1-\frac{N}{K}\right)$$

 - population $N$ grows (or shrinks) with *constant rate* $r$ 
 - only a finite maximum population (*carrying capacity* $K$) can be sustained: limited by resources or other forms of competition among population members
 
The analytic solution is the logistic function with initial population size $N(t=0)=N_0$:

$$N(t) = \frac{N_0\, K}{(K-N_0)\exp(-r\,t)+N_0}=\frac{N_0\, K\exp(r\,t)}{K+ N_0\,(\exp(r\,t)-1)}$$

## Numeric Integration

### Define computational function for ODE

Define a python function that implements the ODE: 

  $$f(t,y)= f(y) = \frac{d\, y}{d\, t}=r\, y \left(1-\frac{y}{K}\right)$$
  
The function should accept input *arguments* `t`, `y`. 
For now, 'hard-code' parameters `r`, `K` in the function body.
 


```python
def ode_logistic_growth(t, y):
    K = 1
    r = 1
    return r*y-r*y*y/K
```

### Solver


```python
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: ODE-defining function object, 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(t.shape)
    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
```

### Solve ODE

- Define array `t` that discretizes time interval $[t_0, t_{\text{max}}]$ in $n$ subintervals of length $\Delta t$.
- Apply numerical ODE solver to your ODE-defining function, providing
    - initial value $N_0$ 
    - array of integration time points $t_i$: `t`
    - other parameters of your ODE as keyword arguments
  Use: $t_0=0\, h$, $t_{\text{max}}=120\, h$, $\Delta t = 0.1$, $N_0=10^3$, $r=0.12 / h$, $K=10^6$
    
- Plot solution.

- Check if initial condition is fulfilled and if solution shows the expected limit behavior.


```python
# ODE with specific parameters
def ode_logistic_growth_with_params(t, y):
    K = 1E6
    r = 0.12
    return r*y-r*y*y/K

# time discretization
t = np.arange(0, 120, 0.1)
# solve ODE
y_numeric = solve_euler(ode_logistic_growth_with_params, t, y_0=1E3)
```


```python
plt.plot(t, y_numeric, label="numeric solution")
plt.legend()
```

## Comparison to Analytic Solution

Testing a new numeric method or implementation on simple test problems with known solutions is a good way to identify problems and gain confidence in the results it produces. 

We will compare the numerical solution of the ODE to its analytical solution.

### Compute analytic solution

Write a function that implements the analytic solution of the logistic growth problem: 


$$N(t) = \frac{N_0\, K}{(K-N_0)\exp(-r\,t)+N_0}=\frac{N_0\, K\exp(r\,t)}{K+ N_0\,(\exp(r\,t)-1)}\, .$$


```python
def logistic_growth(N0, t, r, K):
    enum = N0*K*np.exp(r*t)
    denom = K + N0*np.exp(r*t) - N0
    return enum/denom
```


```python
y_analytic = logistic_growth(1E3, t, r=0.12, K=1E6)
plt.plot(t, y_analytic, label="analytic solution")
plt.legend()
```

### Compare analytic and numeric solution

- Qualitatively compare (plot !) the numeric solutions to the analytic solution for different integration step sizes:
     1. $\Delta t = 0.1$
     2. $\Delta t = 1.0$
     3. $\Delta t = 5.0$
 
- Repeat this for $r=0.5$ 


```python
# ODE with specific parameters
def ode_logistic_growth_with_params(t, y):
    K = 1E6
    r = 0.12
    return r*y-r*y*y/K

for delta_t in [5.0, 1.0, 0.1]:
    t = np.arange(0, 120, delta_t)
    y_numeric = solve_euler(ode_logistic_growth_with_params, t, y_0=1E3)
    plt.plot(t, y_numeric, label="numeric solution: $\Delta t = %.2f$"%delta_t)
    
t_analytic = np.arange(0, 120, 0.1)
y_analytic = logistic_growth(1E3, t, 0.12, 1E6)
plt.plot(t, y_analytic, label="analytic solution")
plt.legend()
```

## Logistic Growth with non-constant Growth Rate

So far we have assumed that the growth rate $r$ is constant in time, but more generally $r=r(t)$ may be a function of time. 

If an analytic expression is known for the functional form of the time dependency, the time-dependency can be included directly into the ODE-defining function and the time-dependent problem can be solved in the same way as the time-independent problem.

Consider the following example of a time dependent growth rate (from lecture notes, section 3.4):

\begin{equation}
    r(t)=\frac{r_{in}+r_{fin}}{2} + \frac{r_{fin}-r_{in}}{2}\tanh(t-D),
\end{equation}

with parameters $r_{in}$, $r_{fin}$, and $D$.


```python
def ode_log_growth_time_dep(t, y):
    rin  = -1
    rfin = 1
    K = 3
    D = 2
    r = 0.5*(rin+rfin)+0.5*(rfin-rin)*np.tanh(t-D)
    return r*y-r*y*y/K

t = np.arange(0, 20, 0.1)
y_numeric = solve_euler(ode_log_growth_time_dep, t, y_0=0.5)

plt.plot(t, y_numeric, label="numeric solution")
plt.legend()
```

# Standard Solver Interface

The `scipy.integrate` package provides multiple ODE solvers, most of them can be accessed via the 
[`scipy.integrate.solve_ivp`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp) interface.
By default, it uses the [`scipy.integrate.RK45`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.RK45.html#scipy.integrate.RK45) solver, an explicit *adaptive* Runge-Kutta method.

When using our own solver implementation we could decide how to call the ODE-defining function from within the solver, i.e. which arguments the ODE-defining function accepts and their order. 
Now, to work with `solve_ivp`, the ODE-defining function *must* have the signature `f(t, y)` where `t` is the current evaluation point and `y` the current function value.

The `solve_ivp` interface returns not only an array of estimated function values (as our implementation of the euler method) but instead a results object that contains various types of information. 
You can inspect this object by printing or using `dir()`.


```python
from scipy.integrate import solve_ivp

def ode_logistic_growth_with_params(t, y):
    K = 1E6
    r = 0.12
    return r*y-r*y*y/K

t_0 = 0
t_N = 120
N0 = 1E3

sol = solve_ivp(ode_logistic_growth_with_params, [t_0, t_N], [N0])   
#print("results object: ",sol)
```


```python
plt.plot(sol.t, sol.y[0]) # note that the result y is returned as nested array
```

## Adaptive Solver Settings

`solve_ivp` uses an *adaptive solver*, i.e. the stepsize $\Delta t$ is controlled by estimating the local error and a target accuracy. See 'Adaptive Stepsize' in the [ODE notebook](https://github.com/cohmathonc/biosci670/blob/master/IntroductionComputationalMethods/05_IntroCompMethods_SolvingODEs.ipynb).

- target accuracy can be set by parameters `rtol`, `atol`.
- no limit on maximum stepsize by default, can be imposed via `max_step`


```python
# specific ODE function version for selected parameter sets
def ode_logistic_growth_with_params(t, y):
    K = 1E6
    r = 0.12
    return r*y-r*y*y/K

sol = solve_ivp(fun=ode_logistic_growth_with_params, 
                t_span=[t_0, t_N], y0=[N0], 
                #max_step=5,
                rtol=1E-6
                )   
plt.plot(sol.t, sol.y[0]) 
plt.plot(sol.t, sol.y[0], 'xr') 
```

## Function Evaluation at Specific (Time-)Points

The `t_eval` argument allows to obtain integration results at specific evaluation points.

In contrast to the constant-stepsize solver that we implemented ourselves, `solve_ivp` will integrate the ODE in the entire domain and according to the specified accuracy.
It will only return those function values that correspond to the evaluation points specified by `t_eval`.


```python
sol_selected_points = solve_ivp(fun=ode_logistic_growth_with_params, 
                                 t_span=[t_0, t_N], y0=[N0], 
                                 t_eval=[0, 20, 50, 62.34567, 70, 100],
                                 rtol=1E-3)
plt.plot(sol.t, sol.y[0]) 
plt.plot(sol_selected_points.t, sol_selected_points.y[0], 'sg')
```

# Optional function arguments: Parameters `K`, `r`

We are often interested in investigating the solution of an ODE in function of model parameters, such as `K` and `r` in the case of the logistic growth model.

Until now, we 'hardcoded' these parameters in the definition of the ODE function, and have redefined this function (of defined a new function with distinct name) for each model parametrization:


```python
#  'Hard-code' parameters in the ODE-defining function
def ode_logistic_growth_hardcoded_params(t, y):
    K = 1E6
    r = 0.12
    return r*y-r*y*y/K

sol = solve_ivp(fun=ode_logistic_growth_hardcoded_params, 
                t_span=[t_0, t_N], y0=[N0])   
```

This approach is straightforward for a small number of parameter sets but becomes inconvenient if more than a few parameter combinations have to be explored manually and is not suited for automated exploration of the parameter space (e.g. for data fitting). 

In the following sections we will explore a few alternative approaches that enable us to pass parameters through the ODE solver to the ODE-defining function.
For all of those approaches, we first need a more general ODE function that accepts the parameters of interest as arguments to its function call.


```python
# a general ODE-defining function with parameters as arguments
def ode_logistic_growth_args(t, y, r, K):
    return r*y-r*y*y/K
```

As we know, the ODE solver makes repeated calls to the ODE function during the iterative solution process.
Usually, the ODE solver passes the current time `t` and function value `y` to the ODE function.

We now need to find a way to pass these additional parameters (`r`, `K`) to the ODE function whenever it is called by the ODE solver.

## Case 1: Your own solver -- full control over expected function signature

If you use a solver that you implemented yourself, you can control how the ODE function is called from within the solver.
This allows you to adapt your solver in such a way that it passes any specified model parameters to the ODE function.

### Passing each parameter separately through solver function to ODE

For example you could write a specific version of the solver that can solve logistic growth problems with parameters `r` and `K` by adding these parameters as arguments in the solver function definition (`param_r`, `param_K`), and passing these argument to the ODE-defining function when it is called.


```python
def ode_log_growth(t, y, r, K):
    return r*y-r*y*y/K

def solve_euler_with_params(f, t, y_0, param_r, param_K):
    y = np.zeros(t.shape)
    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], param_r, param_K)
    return y

t = np.arange(0,100,1)
y = solve_euler_with_params(ode_log_growth, t, 0.1, 0.12, 1E6)
```

This approach implies that a new dedicated solver function needs to be defined for ODEs with different numbers of parameters, i.e. a `solve_euler_1param(f, t, y_0, param1)` for ODEs with one parameter, a `solve_euler_2param(f, t, y_0, param1, param2)` for ODEs with two parameters,  a `solve_euler_3param(f, t, y_0, param1, param2, param3)` for ODEs with three parameters, etc.

### Passing an arbitrary number of parameters through solver function to ODE

Instead of passing single parameters separately, we could collect parameters in a *container*, e.g. a `dict` or `list` object, and pass this object to the ODE function.


```python
# 1) parameters as list

def ode_log_growth(t, y, param_list):
    r = param_list[0]
    K = param_list[1]
    return r*y-r*y*y/K

def solve_euler_with_param_list(f, t, y_0, parameters):
    y = np.zeros(t.shape)
    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], parameters)
    return y

t = np.arange(0,100,1)

parameter_list = [0.12, 1E6]
y = solve_euler_with_param_list(ode_log_growth, t, 0.1, parameter_list)
```


```python
# 2) parameters as dictionary

def ode_log_growth(t, y, param_dict):
    r = param_dict['r']
    K = param_dict['K']
    return r*y-r*y*y/K

def solve_euler_with_param_dict(f, t, y_0, parameters):
    y = np.zeros(t.shape)
    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], parameters)
    return y

t = np.arange(0,100,1)

parameters_dict = {'r':0.12, 'K':1E6}
y = solve_euler_with_param_dict(ode_log_growth, t, 0.1, parameters_dict)
```

Note that in both cases (list, dict), the solver function is identical and can handle ODE-defining functions with any arbitrary number of parameters.

However, you must be sure that parameters passed through the solver function are correctly interpreted by the respective ODE-defining function:
- If you have designed the ODE function to expect a parameter *dictionary* with certain keys, then the parameter object passed to the solver function must be of type dictionary and must contain *the keys expected* by the ODE function.
- Likewise, if you have designed the ODE function to expect a parameter *list* of certain length, then the parameter object passed to the solver function must be of type list and its content must be *ordered as expected* by the ODE function.

## Case 2: No control over function signatures


All of the above approaches required us to change the signature of the ODE function call in the definition of the solver function.  

When using an 'external' function, such as `solve_ivp`, this is not easily possible.

Very often, functions that take other functions as arguments and call these functions internally are designed such that they allow additional parameters to be passed into this internal function call. 
This is often implemented by exposing an (optional) argument that accepts a list of parameter values, exactly as we have discussed above in case 1.

Unfortunately, `solve_ivp` does not provide such a functionality; internally it calls the ODE-defining function with signature `f(t, y)` without the possibility to pass any other arguments to `f`.

The workaround in those cases exploits so-called [*lambda expressions*](https://docs.python.org/3/tutorial/controlflow.html#lambda-expressions), which provide a method to define anonymous functions in python.
( This is similar to the [@-handle](https://www.mathworks.com/help/matlab/matlab_prog/creating-a-function-handle.html) in Matlab)

Lambda expressions can be used for inline definition of simple functions (limited to a single expression); consider the following example:


```python
# 'lambda expressions' allow definition of single-expression functions
a=lambda x: 1+x
print(a(1))

#or even:
(lambda x: 1+x)(2)
```

We can exploit this mechanism to replace the ODE-defining function (with arguments `t`, `y`, and other paramers) by an anonymous 'wrapper' function (with arguments `t`, `y`) directly in the `solve_ivp` function call:


```python
def ode_log_growth(t, y, r, K):
    return r*y-r*y*y/K

# using a lambda expression, we can redefine the ODE function directly in the solve_ip function call
sol = solve_ivp(fun=lambda t, y: ode_log_growth(t, y, 0.12, 1E6), 
                t_span=[0, 120], y0=[1E3], rtol=1E-6)   

plt.plot(sol.t, sol.y[0]) 
```

This is the most flexible approach. It can be applied to any of the solvers that we used and for any model parameter specification method.

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

331e09ab5dcc51a4a1176dab383b9d4a675b6b38
ipynb
Jupyter Notebook
GrowthModels/GrowthModels_LogGrowthNumerical.ipynb
sbranciamore/biosci670coh
92e91aa1e0266b54b0831537b857796eaf0b949c
CC-BY-4.0
2019-02-19T21:20:42.000Z
2021-01-12T20:25:42.000Z
GrowthModels/GrowthModels_LogGrowthNumerical.ipynb
danielabler/biosci670
a9fc3e3062cc7ee4f2da3c26a4c5fc22c5cb8cb4
CC-BY-4.0
GrowthModels/GrowthModels_LogGrowthNumerical.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