```python
%matplotlib inline

import matplotlib.pyplot as plt

import math
import numpy as np

from scipy.integrate import solve_ivp
```

- https://docs.scipy.org/doc/scipy-1.3.0/reference/integrate.html#solving-initial-value-problems-for-ode-systems
- https://docs.scipy.org/doc/scipy-1.3.0/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp (the new API)

- https://docs.scipy.org/doc/scipy-1.3.0/reference/tutorial/integrate.html#ordinary-differential-equations-odeint (the old API)

## Conversion d'une ODE d'ordre $n$ en un système d'ODE d'ordre 1

C.f. "Toutes les mathématiques et les bases de l'informatique" de H.Stocker (Dunod) p.745

### Exemple 1

$$y'' + 5y'^2 + 7y^5 = 6x^3$$

On pose $y_1 = y'$ et donc $y_1' = y''$.

On obtient alors le système couplé:

$$
\left\{
    \begin{array}{ll}
        y_1' + 5 y_1^2 + 7 y^5 & = 6 x^3 \\
        y' - y_1 & = 0
    \end{array}
\right.
$$



### Exemple 2

$$\color{red}{y^{(6)}} + 8 x^2 \color{green}{y^{(4)}} + \ln y = 0$$

On pose :
- $y_1 = y'$
- $y_2 = y_1' = y''$
- $y_3 = y_2' = y_1'' = y'''$
- $\color{green}{y_4} = y_3' = \dots = \color{green}{y^{(4)}}$
- $y_5 = y_4' = \dots = y^{(5)}$
- $\color{red}{y_6} = y_5' = \dots = \color{red}{y^{(6)}}$

On obtient alors le système couplé:

$$
\left\{
    \begin{array}{ll}
        \color{red}{y_6} + 8x^2 \color{green}{y_4} + \ln y & = 0 \\
        y_1 - y' & = 0 \\
        y_2 - y_1' & = 0 \\
        y_3 - y_2' & = 0 \\
        y_4 - y_3' & = 0 \\
        y_5 - y_4' & = 0 \\
        y_6 - y_5' & = 0 \\
    \end{array}
\right.
$$


## Solving initial value problems (IVP) for ODE systems

- https://pundit.pratt.duke.edu/wiki/Python:Ordinary_Differential_Equations
- https://pundit.pratt.duke.edu/wiki/Python:Ordinary_Differential_Equations/Examples

### Example 1

Taken from https://pundit.pratt.duke.edu/piki/index.php?title=Python:Ordinary_Differential_Equations&printable=yes#Simple_Example

Solve

$$y' + y = t$$

i.e.

$$\frac{\text{d}}{\text{d}t} y(t) = t - y(t)$$

for some times between 0 and 15, assuming the initial value for $y(0)$ is 2.


```python
integration_range = [0, 15]
initial_values = [2]

# Define the dependent function f
def f(t, y):
    dydt = t-y
    return dydt

sol = solve_ivp(f, integration_range, initial_values)
print(sol)

plt.plot(sol.t, sol.y[0], 'o-');
```


```python
integration_range = [0, 15]
initial_values = [2]

# Define the dependent function f
def f(t, y):
    dydt = t-y
    return dydt

sol = solve_ivp(f, integration_range, initial_values, dense_output=True)
print(sol)

t = np.linspace(0, 15, 50)
y = sol.sol(t)               # <- with dense_output=True, solve_ivp() returned a callable (function) solution instead of a sample (array)

plt.plot(t, y[0], 'o-');
```


```python
integration_range = [0, 15]
initial_values = [2]

# Define the dependent function f
def f(t, y):
    dydt = t-y
    #print(y.shape)
    return dydt

sol = solve_ivp(f, integration_range, initial_values,
                t_eval = np.linspace(0, 15, 30))         # <--
print(sol)

plt.plot(sol.t, sol.y[0], 'o-');
```

### Example 2

Taken from https://pundit.pratt.duke.edu/piki/index.php?title=Python:Ordinary_Differential_Equations&printable=yes#Another_Example

Solve

$$\frac{\text{d}}{\text{d}t} y(t) + \frac{1}{4} y(t) - x(t) = 0$$

i.e.

$$\frac{\text{d}}{\text{d}t} y(t) = - \frac{1}{4} y(t) + x(t)$$

assuming

$$
\begin{align}
x(t) & = \cos(3t) \\
y(0) & = -1
\end{align}
$$

for some times between 0 and 15.


```python
integration_range = [0, 15]
initial_values = [-1]

def f(t, y):
    dydt = -y + np.cos(3. * t)
    return dydt

sol = solve_ivp(f, integration_range, initial_values,
                t_eval = np.linspace(0, 15, 100))
print(sol)

plt.plot(sol.t, sol.y[0], 'o-');
```

### Example 3 : constant rate of change

Taken from https://pundit.pratt.duke.edu/piki/index.php?title=Python:Ordinary_Differential_Equations/Examples&printable=yes#Constant_Rate_of_Change

Solve

$$\frac{\text{d}}{\text{d}t} y(t) = c$$

assuming

$$
y(0) = 6
$$

for some times between 0 and 10.


```python
integration_range = [0, 10]
initial_values = [6]
c = 1.2

# Define the dependent function f
def f(t, y):
    dydt = [c]
    return dydt

sol = solve_ivp(f, integration_range, initial_values)
print(sol)

plt.plot(sol.t, sol.y[0], 'o-');
```

### Example 4 : time-dependent rate of change

Taken from https://pundit.pratt.duke.edu/piki/index.php?title=Python:Ordinary_Differential_Equations/Examples&printable=yes#Time-dependent_Rate_of_Change

Solve

$$\frac{\text{d}}{\text{d}t} y(t) = 2 t^2 - 6 t + 3$$

assuming

$$
y(0) = 6
$$

for some times between 0 and 4.


```python
integration_range = [0, 4]
initial_values = [6]

# Define the dependent function f
def f(t, y):
    dydt = 2. * t**2 - 6. * t + 3.
    return dydt

sol = solve_ivp(f, integration_range, initial_values,
                t_eval=np.linspace(0, 4, 30))
print(sol)

plt.plot(sol.t, sol.y[0], 'o-');
```

### Example 5 : population growth

Taken from https://pundit.pratt.duke.edu/piki/index.php?title=Python:Ordinary_Differential_Equations/Examples&printable=yes#Population_Growth

For population growth, the rate of change of population is dependent upon the number of people as well as some constant of proportionality:

$$\frac{\text{d}}{\text{d}t} y(t) = c y(t)$$

assuming

$$
\begin{align}
y(0) & = 10   \\
c    & = 1.02
\end{align}
$$

for some times between 0 and 3.


```python
integration_range = [0, 3]
initial_values = [10]
c = 1.02

# Define the dependent function f
def f(t, y):
    dydt = c * y
    return dydt

sol = solve_ivp(f, integration_range, initial_values,
                t_eval=np.linspace(0, 3, 30))
print(sol)

plt.plot(sol.t, sol.y[0], 'o-');
```

### Example 6 : multiple dependent variable models

Taken from https://pundit.pratt.duke.edu/piki/index.php?title=Python:Ordinary_Differential_Equations/Examples&printable=yes#Multiple_Variable_Models

Solve:

$$
\left\{
\begin{align}
\frac{\text{d}}{\text{d}t} y_0(t) & = 4 \cos(3t) \\
\frac{\text{d}}{\text{d}t} y_1(t) & = -2 y_0(t) + 0.5 t
\end{align}
\right.
$$

assuming

$$
\begin{align}
y_1(0) & = 0   \\
y_2(0) & = -3
\end{align}
$$

for some times between 0 and 5.


```python
integration_range = [0, 5]
initial_values = [0, -3]

# Define the dependent function f
def f(t, y):
    dy1dt = 4. * math.cos(3. * t)
    dy2dt = -2. * y[0] + 0.5 * t
    return [dy1dt, dy2dt]

sol = solve_ivp(f, integration_range, initial_values,
                t_eval=np.linspace(0, 5, 50))
print(sol)

fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(16, 4))
ax1.plot(sol.t, sol.y[0], 'o-', label="y1")
ax1.legend()
ax2.plot(sol.t, sol.y[1], 'o-', label="y2")
ax2.legend();
```

### Example 7 : higher Order Differential Equations

Taken from https://pundit.pratt.duke.edu/piki/index.php?title=Python:Ordinary_Differential_Equations/Examples&printable=yes#Higher_Order_Differential_Equations

Solve:

$$
\frac{\text{d}^3}{\text{d}t^3} y(t) = c
$$

This third-order equation can be rewritten as a system of three first-order differential equations:

$$
\left\{
\begin{align}
y_0 = y \\
\frac{\text{d}}{\text{d}t}y_0(t) & = \frac{\text{d}}{\text{d}t}y(t) = y_1 \\
\frac{\text{d}}{\text{d}t}y_1(t) & = \frac{\text{d}^2}{\text{d}t^2} y_0(t) = \frac{\text{d}^2}{\text{d}t^2} y(t) = y_2 \\
\frac{\text{d}}{\text{d}t}y_2(t) & = \frac{\text{d}^2}{\text{d}t^2} y_1(t) = \frac{\text{d}^3}{\text{d}t^3} y_0(t) = \frac{\text{d}^3}{\text{d}t^3} y(t) = c \\
\end{align}
\right.
$$
 
$$
\Leftrightarrow
\left\{
\begin{align}
y_0 = y \\
\frac{\text{d}}{\text{d}t} y_0(t) & = y_1 \\
\frac{\text{d}}{\text{d}t} y_1(t) & = y_2 \\
\frac{\text{d}}{\text{d}t} y_2(t) & = c \\
\end{align}
\right.
$$

assuming

$$
\begin{align}
y_1(0) & = 6   \\
y_2(0) & = 2   \\
y_3(0) & = -4
\end{align}
$$

i.e.

- an initial position of 6
- an initial velocity of 2
- an initial acceleration of -4

and a constant acceleration of 1.3
over a period of 8 seconds.


```python
integration_range = [0, 8]
initial_values = [6, 2, -4]
c = 1.3

# Define the dependent function f
def f(t, y):
    dy1dt = y[1]
    dy2dt = y[2]
    dy3dt = c
    return [dy1dt, dy2dt, dy3dt]

sol = solve_ivp(f, integration_range, initial_values,
                t_eval=np.linspace(0, 8, 50))
print(sol)

fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(14, 4))
ax1.plot(sol.t, sol.y[0], 'o-', label="y1")
ax1.legend()
ax2.plot(sol.t, sol.y[1], 'o-', label="y2")
ax2.legend();
ax3.plot(sol.t, sol.y[2], 'o-', label="y3")
ax3.legend();
```

### Remark about the initial value

Taken from https://pundit.pratt.duke.edu/piki/index.php?title=Python:Ordinary_Differential_Equations&printable=yes#Notes_.2F_Troubleshooting

Lets imagine we want to solve

$$\frac{\text{d}}{\text{d}t} y(t) = c$$

assuming

$$
y(0) = 1
$$

for some times between **5** and 10.


```python
integration_range = [5, 10]
initial_values = [1]
c = 1.2

# Define the dependent function f
def f(t, y):
    dydt = [c]
    return dydt

res = solve_ivp(f, integration_range, initial_values, t_eval=np.linspace(5, 10, 10))
print(res)

plt.plot(res.t, res.y[0], 'o-');
```

Here, we can see that scipy considers that the actual initial value $y(0) = 1$ is not for $t=0$ but for the start of the time span !
