# ODE Explorer Demo 1: Exponential decay models

A small appetizer example to find your way around this library. Introduces the most central components. Let's get to it!


```python
import ode_explorer
from ode_explorer.models import ODEModel
from ode_explorer.integrators import Integrator
from ode_explorer.stepfunctions import RungeKutta4, ForwardEulerMethod, DOPRI45
from ode_explorer.stepsize_control import DOPRI45Controller
from ode_explorer.metrics import DistanceToSolution

import numpy as np

from typing import Union
```

## Model definition

We start with the simplest of models - exponential decay. But do not assume that just because it is simple, it is just a toy example! In fact, the exponential decay model plays an important role in stability analysis of numerical integration methods, where it is also known as the *Dahlquist test equation*.

The differential equation for exponential decay with a rate $\lambda$ reads:
\begin{align}
y' &= - \lambda \cdot y, \\
y(0) &= y_0.
\end{align}

This we can easily cast into a model:


```python
def exp_decay(t: float, y: Union[float, np.ndarray], lamb: float):
    return - lamb * y
```

Small side note: Due to clever operator overloading by NumPy, you can use the above model both for a scalar ODE, i.e. $y\in\mathbb{R}$, as well as for any higher dimension $y\in\mathbb{R}^n$. In fact, you can do the same with the decay parameter $\lambda$ - you can initialize it as a NumPy array $\lambda\in\mathbb{R}^n$ to simulate exponential decay with many different growth rates at the same time!

## Initial state

For any initial value problem (IVP), you need to provide a starting value of the process - that's why it is named that way! Try out different values by editing the cells below.


```python
t_0 = 0.0
y_0 = 1.0  # np.ones(100)
lamb = 0.5 # decay rate of 0.5, so we expect half-life to be 2 * ln(2) ~ 1.38


initial_state = (t_0, y_0)
```

## Model construction

We wrap our exponential decay model into an ODEModel class.


```python
model = ODEModel(ode_fn=exp_decay, fn_args={"lamb": lamb})
```

## Integrator construction

To solve our exponential decay differential equation, we instantiate an *Integrator* object, which keeps some useful state around for us.


```python
integrator = Integrator()
```

    I1230 19:44:13.571961 4548111872 integrator.py:98] Creating an Integrator instance.
    I1230 19:44:13.573229 4548111872 integrator.py:83] Created an Integrator instance.


## Integration

That's how easy it is. Now let's solve our equation!


```python
integrator.integrate_const(model=model,
                           step_func=RungeKutta4(),
                           initial_state=initial_state,
                           h=0.001,
                           max_steps=10000,
                           verbosity=1,
                           progress_bar=True)
```

    I1230 19:44:13.579713 4548111872 integrator.py:191] Starting integration.
    100%|██████████| 10000/10000 [00:00<00:00, 62566.44it/s]
    I1230 19:44:13.747277 4548111872 integrator.py:200] Finished integration.





    <ode_explorer.integrators.integrator.Integrator at 0x12e1b1c90>



As we can see from the progress bar, in the scalar case we managed about 24,000 iterations per second with the classical fourth order Runge-Kutta method. A successful numerical integration is called a *run* - let's find our experiment first by listing all runs.


```python
integrator.list_runs()
```

    | timestamp                | run_id                               |
    |--------------------------|--------------------------------------|
    | Wed Dec 30 19:44:13 2020 | 78bd4868-fffd-4894-bc39-283601669562 |


And there we have our runs, neatly tabulated. The more experiments you do, the more runs will show up in this table. Each run has its own identifier called *run_id*, which is a UUID generated during the run. 

We also see that the output format of our step function is "zipped" - that means that internally, the step function returns a dictionary of floats instead of a dictionary of state vectors. This is beneficial for later analysis with pandas, but creates a performance hit. Try out the other option if you want - just supply output_format="variables" to the RungeKutta4 constructor.

Let us look at the result of our first run - instead of selecting it by the run's ID as you normally would, we can also type "latest" to indicate that we want the newest result.


```python
result = integrator.return_result_data(run_id="latest")

result
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>t</th>
      <th>y</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0.000</td>
      <td>1.000000</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0.001</td>
      <td>0.999500</td>
    </tr>
    <tr>
      <th>2</th>
      <td>0.002</td>
      <td>0.999000</td>
    </tr>
    <tr>
      <th>3</th>
      <td>0.003</td>
      <td>0.998501</td>
    </tr>
    <tr>
      <th>4</th>
      <td>0.004</td>
      <td>0.998002</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
      <th>9996</th>
      <td>9.996</td>
      <td>0.006751</td>
    </tr>
    <tr>
      <th>9997</th>
      <td>9.997</td>
      <td>0.006748</td>
    </tr>
    <tr>
      <th>9998</th>
      <td>9.998</td>
      <td>0.006745</td>
    </tr>
    <tr>
      <th>9999</th>
      <td>9.999</td>
      <td>0.006741</td>
    </tr>
    <tr>
      <th>10000</th>
      <td>10.000</td>
      <td>0.006738</td>
    </tr>
  </tbody>
</table>
<p>10001 rows × 2 columns</p>
</div>



That looks ok - although staring at the numbers is not really the best option here for assessing whether our solution actually makes sense. If you have worked with pandas before, you may have noticed by the print output that the returned "result" object is in fact a pandas DataFrame - which can be plotted immediately. Let's do it!


```python
result.plot(x="t", y="y")
```

That looks much more familiar! Great job, the general shape of the curve looks correct for our example. But how accurate is the solution really? That is something we cannot directly infer from this graph. And if you had to assert numbers like these for a Mars mission launch, you would surely want a way to validate that your algorithm does not go crazy - a lot could be at stake!

Luckily, for an exponential decay model like this, we have a closed form solution available - literally the exponential function. So how about we check our calculations against the real solution?

## Solution distance experiment

We define our solution, which mathematically is a function $y: \mathbb{R}\rightarrow\mathbb{R}^n$. For the exponential decay model with decay rate $\lambda$ and initial state $y_0\in\mathbb{R}^n$, the exact solution is given for every $t \geq 0$ as
\begin{equation}
\hat{y}(t) = y_0 \exp(-\lambda t).
\end{equation}


```python
def sol(t):
    return y_0 * np.exp(-lamb * t)
```

Nice! Notice we reused the $y_0$ object from earlier in the notebook. Now we start another integrator run with the same parameters, but this time we add a *Metric*, the *Euclidean distance from the exact solution* $\hat{y}$
\begin{equation}
d(y, \hat{y}, t) = \left\lVert y(t) - \hat{y}(t) \right\rVert_2 = \sqrt{\sum_i (y_i - \hat{y}_i)^2}.
\end{equation}

If our algorithm is correct, which it should be (it is a built-in!), the distance should be small throughout the run, indicating an accurate numerical approximation. Fingers crossed!

*Tip: You can also calculate other norms, like the maximum norm or the $L_1$ distance between solution and approximation. To try out other norms, pass the **norm** argument in the DistanceToSolution constructor (norm=1 is $L_1$ distance, norm="inf" is maximum norm distance)*


```python
integrator.integrate_const(model=model,
                           step_func=RungeKutta4(),
                           initial_state=initial_state,
                           h=0.001,
                           max_steps=10000,
                           verbosity=1,
                           progress_bar=True,
                           metrics=[DistanceToSolution(solution=sol, name="l2_distance")])
```

    I1230 19:44:14.256709 4548111872 integrator.py:191] Starting integration.
    100%|██████████| 10000/10000 [00:00<00:00, 31518.39it/s]
    I1230 19:44:14.577193 4548111872 integrator.py:200] Finished integration.





    <ode_explorer.integrators.integrator.Integrator at 0x12e1b1c90>



Looks like that worked! Let us quickly confirm.


```python
integrator.list_runs()
```

    | timestamp                | run_id                               |
    |--------------------------|--------------------------------------|
    | Wed Dec 30 19:44:13 2020 | 78bd4868-fffd-4894-bc39-283601669562 |
    | Wed Dec 30 19:44:14 2020 | 97466b81-d7ab-4695-b00f-860fb9d6ff9f |


Perfect! There is a new run listed in the table. It should be holding the distance to the solution as metric values. To request the metric data, we can use the *return_metrics* interface exposed by the integrator.


```python
metrics = integrator.return_metrics(run_id="latest")

metrics
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>l2_distance</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0.000000e+00</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0.000000e+00</td>
    </tr>
    <tr>
      <th>2</th>
      <td>0.000000e+00</td>
    </tr>
    <tr>
      <th>3</th>
      <td>0.000000e+00</td>
    </tr>
    <tr>
      <th>4</th>
      <td>0.000000e+00</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
    </tr>
    <tr>
      <th>9996</th>
      <td>3.417405e-16</td>
    </tr>
    <tr>
      <th>9997</th>
      <td>3.426079e-16</td>
    </tr>
    <tr>
      <th>9998</th>
      <td>3.443426e-16</td>
    </tr>
    <tr>
      <th>9999</th>
      <td>3.460773e-16</td>
    </tr>
    <tr>
      <th>10000</th>
      <td>3.478121e-16</td>
    </tr>
  </tbody>
</table>
<p>10001 rows × 1 columns</p>
</div>



Lots of zeros and numbers around machine precision there - that looks promising! The metrics object above is also a pandas DataFrame, so let us plot the solution distance over time:


```python
metrics.plot()
```

From this plot, it appears that the Euclidean distance of our numerical solution to the exact solution never exceeds $3\cdot 10^{-14}$ across all 10,000 iterations. Though the distance does not appear to be uniformly at around machine precision ($\approx 10^{-16}$), but rather periodically increasing and decreasing - an interesting artifact. We can quickly print some statistics, too:


```python
metrics.describe()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>l2_distance</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>count</th>
      <td>1.000100e+04</td>
    </tr>
    <tr>
      <th>mean</th>
      <td>9.958515e-15</td>
    </tr>
    <tr>
      <th>std</th>
      <td>8.205534e-15</td>
    </tr>
    <tr>
      <th>min</th>
      <td>0.000000e+00</td>
    </tr>
    <tr>
      <th>25%</th>
      <td>2.655862e-15</td>
    </tr>
    <tr>
      <th>50%</th>
      <td>9.103829e-15</td>
    </tr>
    <tr>
      <th>75%</th>
      <td>1.618150e-14</td>
    </tr>
    <tr>
      <th>max</th>
      <td>2.473022e-14</td>
    </tr>
  </tbody>
</table>
</div>



That confirms what we were seeing. In theory, the classical Runge-Kutta method has order of consistency 4, meaning that the discretization error is globally of $\mathcal{O}(h^4)$. For a step size of $h = 0.001$ that we chose, this would place us in the range of $10^{-12}$ relative to $h = 1$ (ignoring constants) - which seems all right given our results, and considering that the condition number of the exponential decay equation is not large for $\lambda = 0.5$, no stiffness effects are coming to haunt us here, either. 

This concludes the first small demo of the ode-explorer library - if you made it this far, thanks for trying it out, and see you in the next one!

13879ccf1d34dba9f2edbba4b654f4b5d1808e06
ipynb
Jupyter Notebook
ode_explorer/examples/exponential_decay/exponential_decay.ipynb
njunge94/ode-explorer
0bcb5d4834f6001b2a3e54bd5e000e86bbedf221
MIT
2020-12-21T20:09:59.000Z
2021-03-03T14:27:45.000Z
ode_explorer/examples/exponential_decay/exponential_decay.ipynb
njunge94/ode-explorer
0bcb5d4834f6001b2a3e54bd5e000e86bbedf221
MIT
ode_explorer/examples/exponential_decay/exponential_decay.ipynb
njunge94/ode-explorer
0bcb5d4834f6001b2a3e54bd5e000e86bbedf221
MIT
Qwen/Qwen-72B
1. YES
2. YES