notebook.community

Edit and run 





In [ ]:



    


%matplotlib inline

import matplotlib.pyplot as plt

import math
import numpy as np

from scipy.integrate import solve_ivp

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

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.





In [ ]:



    


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-');



    











In [ ]:



    


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-');



    











In [ ]:



    


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.





In [ ]:



    


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.





In [ ]:



    


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.





In [ ]:



    


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.





In [ ]:



    


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.





In [ ]:



    


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.





In [ ]:



    


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.





In [ ]:



    


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 !

Content source: jdhp-docs/notebooks

Similar notebooks:

notebook.community | gallery | about