Intégration et dérivation numérique¶

Utilisation des modules scipy et numpy pour l'intégration et la dérivation numérique.

Intégration numérique¶

Intégrale d'une fonction¶

Les fonctions d'intégration disponible avec scipy.integrate sont :

Fonction	Application
`quad(f,a,b)`	Calcul une valeur approchée de $\int_a^bf$ par une méthode optimisée par Python.
`dblquad(f,a,b)` et `tplquad(f,a,b)`	Même application que `quad` mais pour des intégrales double et triple.
`romberg(f,a,b)`	Calcul une valeur approchée de $\int_a^bf$ par la méthode de Romberg.

Soit l'intégrale suivante : $I = \int_0^1\frac{1}{1+x}dx$

from scipy.integrate import quad  # Module d'intégration "quad"
import numpy as np

# Intervalle d'intégration
a = 0
b = 1

# Fonction à intégrer
def function(x):
    return 1/(1+x)


# Calcul de l'intégrale
res, err = quad(function, a, b)

# Affichage du résultat
print("Résultat de l'intégrale :", res)

    Résultat de l'intégrale : 0.6931471805599454

from scipy.integrate import romberg  # Module d'intégration "romberg"
import numpy as np

# Intervalle d'intégration
a = 0.
b = 1.

# Fonction à intégrer
def fun(x): return 1/(1+x)


# Calcul de l'intégrale
res = romberg(fun, a, b)

# Affichage du résultat
print("Résultat de l'intégrale :", res)

    Résultat de l'intégrale : 0.6931471805622968

Intégrale double¶

Soit l'intégrale : $I = \iint_0^1x\,y\,dx\,dy$

from scipy.integrate import dblquad  # Module d'intégration "dblquad"
import numpy as np

# Intervalle d'intégration
a = 0
b = 1

# Fonction à intégrer
def function(y, x):
    return x*y


# Calcul de l'intégrale
res, err = dblquad(function, a, b, lambda x: a, lambda x: b)

# Affichage du résultat
print("Résultat de l'intégrale :", res)

    Résultat de l'intégrale : 0.24999999999999997

Soit l'intégrale : $I = \int_0^2\left[\int_0^1x\,y\,dx\right]\,dy$

from scipy.integrate import dblquad  # Module d'intégration "dblquad"
import numpy as np

# Intervalle d'intégration
ax = 0
bx = 1
ay = 0
by = 2

# Fonction à intégrer
def func(y, x): return x*y


# Calcul de l'intégrale : lambda x : ay permet bien de choisir l'intervalle de y
res, err = dblquad(func, ax, bx, lambda x: ay, lambda x: by)

# Affichage du résultat
print("Résultat de l'intégrale :", res)

    Résultat de l'intégrale : 0.9999999999999999

Intégrale triple¶

Soit l'intégrale : $I = \int_0^3\left[\int_0^2\left[\int_0^1x\,y\,z\,dx\right]\,dy\right]dz$

from scipy.integrate import tplquad  # Module d'intégration "tplquad"
import numpy as np

# Intervalle d'intégration
ax = 0
bx = 1
ay = 0
by = 2
az = 0
bz = 3

# Fonction à intégrer
def func(z, y, x): return x*y*z


# Calcul de l'intégrale
res, err = tplquad(func, ax, bx, lambda x: ay, lambda x: by,
                   lambda x, y: az, lambda x, y: bz)

# Affichage du résultat
print("Résultat de l'intégrale :", res)

    Résultat de l'intégrale : 4.5

Intégration d'une liste de valeurs¶

Les fonctions d'intégration disponible avec scipy.integrate sont :

Fonction	Application
`cumtrapz(y,x,initial=0)`	Calcul l'intégrale cumulée par la méthode des trapèzes.
`trapz(y,x)`	Calcul l'intégrale pour chaque point par la méthode des trapèzes.
`simps(y,x)`	Calcul l'intégrale pour chaque point par la méthode de Simpson.

Il existe dans le module numpy, la fonction trapz(y,x) qui permet d'intégrer une liste de valeurs.

Intégrale cumulée¶

L'intégration cumulée permet d'avoir une liste de valeurs contenant le résultat de l'intégrale à chaque instant.

Pour cet exemple, on récupère la vitesse d'une balle en fonction du temps lors d'une chute libre.
Le fichier est un fichier csv se présentant sous la forme suivante :

    t ;x ;y ;v 
    s ;m ;m ;m / s 
    0;0;-0,004;0
    0,04;0;-0,004;0
    0,08;0;-0,004;0
    0,12;0;-0,004;0,39
    0,16;0;-0,036;0,89
    0,2;0;-0,076;1,29
    0,24;0;-0,139;1,68
    0,28;0,004;-0,21;2,06
    0,32;0;-0,304;2,51
    0,36;0,004;-0,411;2,85
    0,4;0,004;-0,532;3,19
    0,44;0;-0,666;3,64
    0,48;-0,004;-0,823;4,03
    0,52;-0,004;-0,988;4,46
    0,56;0;-1,18;4,9
    0,6;0;-1,38;5,25
    0,64;0;-1,6;5,63
    0,68;0;-1,83;6
    0,72;-0,004;-2,08;6,38
    0,76;0;-2,34;

from scipy.integrate import cumtrapz  # Module d'intégration "cumtrapz"
import numpy as np

# Lecture du fichier
# loadtxt ne reconnaissant pas la virgule comme décimale,
# il faut procéder à un post-traitement
data = np.loadtxt("Chute.csv", dtype=np.str, delimiter=';', skiprows=2)
data = np.char.replace(data, ',', '.')

# Extraction des données utiles
t = np.asarray(data[:, 0], dtype=np.float, order='C')
y = np.asarray(data[:, 2], dtype=np.float, order='C')
# Ici la vitesse est issue de la dérivé de la position
v = np.asarray(data[:-1, 3], dtype=np.float, order='C')


y_calc = cumtrapz(v, t[:-1], initial=0)

# Affichage du résultat
print("Résultat de l'intégrale :", y_calc)

    Résultat de l'intégrale : [ 0.      0.      0.     -0.016  -0.052  -0.1035 -0.1705 -0.253  -0.3535
     -0.4675 -0.595  -0.7405 -0.9015 -1.08   -1.276  -1.486  -1.711  -1.951 ]

import matplotlib.pyplot as plt

# Tracé de la position
plt.plot(t[:-1], v, label="Vitesse")
plt.plot(t, y, label="Position issue du fichier", color="red")
plt.scatter(t[:-1], y_calc, label="Position calculée")
plt.ylabel('Amplitude (m/s) et (m)')
plt.xlabel("Temps (s)")
plt.title("Position et vitesse de la balle")
plt.grid(which="both")
plt.legend()
plt.show()

png

Intégrale complète¶

Intégration à partir de deux listes de valeurs correspondant à deux axes (par exemple une liste pour la mesure et une pour le temps le temps).

from scipy.integrate import trapz  # Module d'intégration "trapz"
import numpy as np

# Création d'une liste de valeur à partir de l'équation 1/(1+x) entre 0 et 1
a = 0
b = 1
points = np.linspace(a, b, 100)

y = 1/(1+points)

# Calcul de l'intégral
res = trapz(y, x=points)

# Affichage du résultat
print("Résultat de l'intégrale :", res)

    Résultat de l'intégrale : 0.6931535573789361

from scipy.integrate import simps  # Module d'intégration "simps"
import numpy as np

# Création d'une liste de valeur à partir de l'équation 1/(1+x) entre 0 et 1
a = 0
b = 1
points = np.linspace(a, b, 100)

y = 1/(1+points)

# Calcul de l'intégral
res = simps(y, x=points)

# Affichage du résultat
print("Résultat de l'intégrale :", res)

    Résultat de l'intégrale : 0.6931472762939747

import numpy as np

# Création d'une liste de valeur à partir de l'équation 1/(1+x) entre 0 et 1
a = 0
b = 1
points = np.linspace(a, b, 100)

y = 1/(1+points)

# Calcul de l'intégral
res = np.trapz(y, x=points)

# Affichage du résultat
print("Résultat de l'intégrale :", res)

    Résultat de l'intégrale : 0.6931535573789361

Info

Si vous n'avez pas de deuxième liste par exemple les instants de mesure, mais uniquement l'intervalle entre deux points, vous pouvez passer le paramètre dx=... à la place de x=....
Ce paramètre fonctionne pour l'ensemble des méthodes d'intégration.

from scipy.integrate import trapz  # Module d'intégration "trapz"
import numpy as np

# Création d'une liste de valeur à partir de l'équation 1/(1+x) entre 0 et 1
a = 0
b = 1
points = np.linspace(a, b, 100)
intervalle = 1/99

y = 1/(1+points)

# Calcul de l'intégral
res = trapz(y, x=points)
res2 = trapz(y, dx=intervalle)

# Affichage du résultat
print("Résultat de l'intégrale par liste :", res)
print("Résultat de l'intégrale par intervalle :", res2)

    Résultat de l'intégrale par liste : 0.6931535573789361
    Résultat de l'intégrale par intervalle : 0.6931535573789361

Dérivation numérique¶

Dérivée d'une fonction¶

On utilisera la fonction derivative(f,x0) du module scipy.misc pour calculer la dérivée de fen x0.

Soit la fonction suivante : $f(x) = \frac{1}{1+x}$

Calcul de la dérivée en un point précis¶

from scipy.misc import derivative  # Fonction de dérivation
import numpy as np

# Fonction à intégrer
def function(x):
    return 1/(1+x)


# Point où l'on veut intégrer
x0 = 1

# Calcul de la dérivée
# dx correspond à l'intervalle entre deux points pour le calcul
res = derivative(function, x0, dx=1e-6)

# Affichage du résultat
print("Résultat de la dérivée :", res)

    Résultat de la dérivée : -0.2500000000349445

Calcul de la dérivée pour une série de points¶

On peut aussi chercher les valeurs de la dérivée pour une série de points.

from scipy.misc import derivative  # Fonction de dérivation
import numpy as np

# Fonction à dérivée
def function(x):
    return 1/(1+x)


# Points où l'on veut la dérivée
x = np.linspace(-10, 10, 100)

# Calcul de la dérivée
res = derivative(function, x, dx=1e-6)


# Affichage du résultat
print("Résultat de la dérivée :", res)

    Résultat de la dérivée : [-1.23456790e-02 -1.29191541e-02 -1.35335356e-02 -1.41928084e-02
     -1.49014552e-02 -1.56645319e-02 -1.64877591e-02 -1.73776288e-02
     -1.83415332e-02 -1.93879186e-02 -2.05264712e-02 -2.17683419e-02
     -2.31264202e-02 -2.46156706e-02 -2.62535459e-02 -2.80605014e-02
     -3.00606365e-02 -3.22825024e-02 -3.47601264e-02 -3.75343232e-02
     -4.06543859e-02 -4.41802913e-02 -4.81856038e-02 -5.27613439e-02
     -5.80212052e-02 -6.41086858e-02 -7.12069805e-02 -7.95529257e-02
     -8.94570148e-02 -1.01332699e-01 -1.15740249e-01 -1.33454065e-01
     -1.55568959e-01 -1.83673469e-01 -2.20143303e-01 -2.68660399e-01
     -3.35180055e-01 -4.29849568e-01 -5.71120564e-01 -7.95471147e-01
     -1.18355271e+00 -1.94425709e+00 -3.76816609e+00 -1.01987513e+01
     -8.10000000e+01 -1.21000000e+02 -1.16539834e+01 -4.08204915e+00
     -2.05860113e+00 -1.23734377e+00 -8.24930561e-01 -5.88967009e-01
     -4.41466601e-01 -3.43160254e-01 -2.74376417e-01 -2.24376731e-01
     -1.86895750e-01 -1.58078096e-01 -1.35445889e-01 -1.17347733e-01
     -1.02648695e-01 -9.05479439e-02 -8.04673196e-02 -7.19809637e-02
     -6.47695958e-02 -5.85900371e-02 -5.32544379e-02 -4.86158303e-02
     -4.45578989e-02 -4.09876172e-02 -3.78298679e-02 -3.50234598e-02
     -3.25181403e-02 -3.02723305e-02 -2.82513887e-02 -2.64262661e-02
     -2.47724579e-02 -2.32691755e-02 -2.18986909e-02 -2.06458109e-02
     -1.94974547e-02 -1.84423106e-02 -1.74705571e-02 -1.65736327e-02
     -1.57440472e-02 -1.49752247e-02 -1.42613727e-02 -1.35973729e-02
     -1.29786893e-02 -1.24012901e-02 -1.18615822e-02 -1.13563550e-02
     -1.08827327e-02 -1.04381332e-02 -1.00202327e-02 -9.62693537e-03
     -9.25634727e-03 -8.90675308e-03 -8.57659650e-03 -8.26446281e-03]

import matplotlib.pyplot as plt

y = function(x)
# Tracé de la fonction et sa dérivée
plt.scatter(x, y, label="f(x)", color="red")
plt.plot(x, res, label="f'(x)")
plt.ylabel('y')
plt.xlabel("x")
plt.title("Fonction et sa dérivée")
plt.grid(which="both")
plt.legend()
plt.show()

png

Dérivée d'une liste de point avec pas fixe¶

Pour dériver une liste de point, on utilise la fonction diff(x) du module numpy qui calcul la différence entre deux points successifs et on divise le tout par le pas.

Pour cet exemple on part d'un fichier csv contenant la position d'une balle au cours du temps lors d'une chute libre.
Le fichier est un fichier csv se présentant sous la forme suivante :

    t ;x ;y ;v 
    s ;m ;m ;m / s 
    0;0;-0,004;0
    0,04;0;-0,004;0
    0,08;0;-0,004;0
    0,12;0;-0,004;0,39
    0,16;0;-0,036;0,89
    0,2;0;-0,076;1,29
    0,24;0;-0,139;1,68
    0,28;0,004;-0,21;2,06
    0,32;0;-0,304;2,51
    0,36;0,004;-0,411;2,85
    0,4;0,004;-0,532;3,19
    0,44;0;-0,666;3,64
    0,48;-0,004;-0,823;4,03
    0,52;-0,004;-0,988;4,46
    0,56;0;-1,18;4,9
    0,6;0;-1,38;5,25
    0,64;0;-1,6;5,63
    0,68;0;-1,83;6
    0,72;-0,004;-2,08;6,38
    0,76;0;-2,34;

Attention

Cette méthode ne marche que si l'on a un pas fixe.

import numpy as np

# Lecture du fichier
# loadtxt ne reconnaissant pas la virgule comme décimale,
# il faut procéder à un post-traitement
data = np.loadtxt("Chute.csv", dtype=np.str, delimiter=';', skiprows=2)
data = np.char.replace(data, ',', '.')

# Extraction des données utiles
t = np.asarray(data[:, 0], dtype=np.float, order='C')
y = np.asarray(data[:, 2], dtype=np.float, order='C')
# Ici la vitesse est issue de la dérivé de la position
v = np.asarray(data[:-1, 3], dtype=np.float, order='C')

pas = 0.04

v_calc = np.diff(y)/pas

# Affichage du résultat
print("Résultat de la dérivée :", v_calc)

    Résultat de la dérivée : [ 0.     0.     0.    -0.8   -1.    -1.575 -1.775 -2.35  -2.675 -3.025
     -3.35  -3.925 -4.125 -4.8   -5.    -5.5   -5.75  -6.25  -6.5  ]

import matplotlib.pyplot as plt

# Tracé de la position
plt.plot(t, y, label="Position", color="red")
plt.plot(t[:-1], v, label="Vitesse issue du fichier")
plt.scatter(t[:-1], v_calc, label="Vitesse calculée")
plt.ylabel('Amplitude (m/s) et (m)')
plt.xlabel("Temps (s)")
plt.title("Position et vitesse de la balle")
plt.grid(which="both")
plt.legend()
plt.show()

png

Dérivée d'une liste de point avec pas variable¶

Dans le cas où le pas n'est pas fixe, il faut écrire une fonction qui calcule la dérivée pour chaque point.

import numpy as np

# Lecture du fichier
# loadtxt ne reconnaissant pas la virgule comme décimale,
# il faut procéder à un post-traitement
data = np.loadtxt("Chute.csv", dtype=np.str, delimiter=';', skiprows=2)
data = np.char.replace(data, ',', '.')

# Extraction des données utiles
t = np.asarray(data[:, 0], dtype=np.float, order='C')
y = np.asarray(data[:, 2], dtype=np.float, order='C')
# Ici la vitesse est issue de la dérivé de la position
v = np.asarray(data[:-1, 3], dtype=np.float, order='C')

# Fonction qui calcul la dérivée de deux listes
def derivee(y, x):
    d = []
    for i in range(0, len(y)-1):
        dy = y[i+1]-y[i]
        dx = x[i+1]-x[i]
        d.append(dy/dx)
    return d


v_calc = derivee(y, t)

# Affichage du résultat
print("Résultat de la dérivée :", v_calc)

    Résultat de la dérivée : [0.0, 0.0, 0.0, -0.7999999999999998, -0.9999999999999998, -1.575000000000001, -1.774999999999998, -2.350000000000001, -2.675000000000001, -3.0249999999999986, -3.350000000000002, -3.925, -4.124999999999997, -4.7999999999999945, -5.000000000000008, -5.5, -5.749999999999995, -6.2500000000000115, -6.499999999999989]

import matplotlib.pyplot as plt

# Tracé de la position
plt.plot(t, y, label="Position", color="red")
plt.plot(t[:-1], v, label="Vitesse issue du fichier")
plt.scatter(t[:-1], v_calc, label="Vitesse calculée")
plt.ylabel('Amplitude (m/s) et (m)')
plt.xlabel("Temps (s)")
plt.title("Position et vitesse de la balle")
plt.grid(which="both")
plt.legend()
plt.show()

png