Python pour le calcul scientifique/Calcul différentiel et intégral

Rappelons que dorénavant les programmes commencent tous par :

#!/usr/bin/python3

import numpy as np
import matplotlib.pyplot as plt

Dérivation numérique

Pour une fonction définie

La fonction scipy.misc.derivative() permet de calculer les dérivées numériques successives d'une fonction, par exemple :

import scipy.misc

def f(x):
    return x*x

X = np.array([0, 1, 2, 3, 4])

print(f(X))
# [ 0  1  4  9 16]

print(scipy.misc.derivative(f, X))
# [0. 2. 4. 6. 8.]

print(scipy.misc.derivative(f, X, n=2)) # dérivée seconde
# [2. 2. 2. 2. 2.]

Le module scipy.optimize fournit la fonction approx_fprime() mais celle-ci ne fonctionne qu'avec une unique valeur de x à la fois.

import numpy as np
from scipy import optimize as opt

def carre(X):
    return X[0]*X[0] + X[1]*X[1]

def gradcarre(X):
    return 2*(X[0] + X[1])

x = np.array([0, 0.5])

print(carre(x)) # 0.25

print(opt.approx_fprime(x, carre, 1e-6))
# [9.99977878e-07 1.00000100e+00]

Pour une liste de valeurs

La fonction np.diff() calcule la différence entre les éléments consécutifs d'un vecteur (ou d'une liste ou d'un n-uplet) : np.diff(M) == M[1:] - M[:-1]. Si M est une matrice, il faut indiquer l'axe en paramètre axis= : 0 (premier indice) pour faire une différence entre les lignes, 1 (deuxième indice) pour faire la différence entre les colonnes… L'axe par défaut est -1.

La fonction np.ediff1d() fait la différence en remettant « en ligne » (comme avec .flat) la matrice.

Par exemple :

M = np.array([[1, 2, 3, 4], [3, 4, 5, 6]])

print(np.diff(M, axis=0))
# [[2 2 2 2]

print(np.diff(M, axis=1))
# [[1 1 1]
#  [1 1 1]]

print(np.ediff1d(M))
# [ 1  1  1 -1  1  1  1]

Si les coordonnées des points d'une courbe sont stockées dans deux vecteurs X et Y, les valeurs étant classées par x croissant, alors on peut calculer les taux de variation en divisant les différences, par exemple :

X = np.array([0, 1, 2, 3])
Y = X*X

print(np.diff(Y)/np.diff(X))
# [1. 3. 5.]

La fonction np.gradient() permet de calculer le gradient d'une matrice en faisant localement l'approximation que la fonction est un polynôme de degré 2, le gradient étant alors obtenu par une combinaison linéaire des valeurs consécutives. Ainsi, il y a autant de valeur de gradient que d'éléments de la matrice. Par exemple :

Y = np.array([0, 1, 2, 3])
print(np.gradient(Y))
# [1. 1. 1. 1.]

Z = np.array([[0, 1, 2, 3], [2, 3, 4, 5]])
print(np.gradient(Z))
# [array([[2., 2., 2., 2.],
#         [2., 2., 2., 2.]]),
#  array([[1., 1., 1., 1.],
#         [1., 1., 1., 1.]])
# ]

Par défaut, la fonction considère que l'on a valeurs Y correspondant à des valeurs X espacées de 1 (donc sur une grille carrée dans le cas d'une matrice). On peut indiquer un autre pas homogène : np.gradient(Z, pas). S'il s'agit d'une grille rectangulaire (le pas n'est pas le même selon la dimension), on indique les deux pas : np.gradient(Z, pas1, pas2).

Si la répartition n'est pas régulière, il faut indiquer les coordonnées dans des matrices.

Intégration numérique

Pour une fonction définie

Le module integrate de Scipy fournit plusieurs fonctions pour faire l'intégration numérique d'une fonction définie :

scipy.integrate.quad() : par la méthode de quadrature, donne la valeur de l'intégrale et une estimation de l'erreur ;
scipy.integrate.fixed_quad() : par une méthode de quadrature de Gauss d'ordre fixe, donne la valeur de l'intégrale et none (valeur fixe) ;
scipy.integrate.quadrature() : par une méthode de quadrature de Gauss adaptative, donne la valeur de l'intégrale et la différence entre les deux dernières aleurs calculées ;
scipy.integrate.romberg() : par la méthode de Romberg, donne la valeur de l'intégrale ;
…

from scipy import integrate as intg

def carre(x):
    return x*x

Y = intg.quad(carre, 0, 4)
Y1 = intg.fixed_quad(carre, 0, 4)
Y2 = intg.quadrature(carre, 0, 4)
Y3 = intg.romberg(carre, 0, 4)

print("quad :", Y, "\n", "fixed_quad :", Y1, "\n", , "quadrature :", Y2, "\n", "romberg :", Y3)
# quad : (21.333333333333336, 2.368475785867001e-13) 
# fixed_quad : (21.33333333333333, None) 
# quadrature : (21.333333333333332, 3.552713678800501e-15) 
# romberg : 21.333333333333332

Notons que la valeur exacte vaut :

\int _{0}^{4}x^{2}\mathrm {d} x={\frac {1}{3}}\left[\mathrm {X} ^{3}\right]_{0}^{4}={\frac {4^{3}-0}{3}}={\frac {64}{3}}=21,3{\underline {3}}\ldots

Pour une liste de valeurs

La fonction numpy.sum(), ainsi que la méthode .sum() pour les matrices, permet d'avoir une estimation de l'intégrale par la méthode du point milieu. Si les valeurs de y correspondent à des valeurs de x équidistantes d'une largeur h, alors une première approche de l'intégrale est :

def carre(x):
    return x*x

x = np.array([0, 1, 2, 3, 4])
h = 1
y = carre(x)

I = h*y.sum()
print(I) # 30

La fonction numpy.trapz() fait l'intégrale par la méthode des trapèzes. Si les points sont espacés régulièrement, on utilise le paramètre dx ; sinon, il faut fournir la liste des valeurs x correspondantes. Avec l'exemple précédent :

I1 = np.trapz(y, dx=h)
I2 = np.trapz(y, x)
print(I1, I2) # 22.0 22.0

Le module scipy.integrate fournit plusieurs méthodes pour affiner le calcul de l'intégrale :

scipy.integrate.trapz() : méthode des trapèzes ;
scipy.integrate.simps() : méthode de Simpson (Newton-Cotes pour 3 points) ;
scipy.integrate.romb() : méthode de Romberg ;
…

from scipy import integrate as intg

def carre(x):
    return x*x

x = np.array([0, 1, 2, 3, 4])
y = carre(x)
h = 1

I11 = intg.trapz(y, dx=h) # méthode des trapèzes
I12 = intg.trapz(y, x)
I21 = intg.simps(y, dx=h) # méthode de Simpson (Newton-Cotes pour 3 points)
I22 = intg.simps(y, x)
I3 = intg.romb(y, dx=h) # méthode de Romberg

print("dx, trapèzes :", I11, "\n", "Simpson :", I21, "\n", "Romberg :", I3)
print("x, trapèzes :", I12, "\n", "Simpson :", I22, "\n")
# dx, trapèzes : 22.0 
#  Simpson : 21.333333333333332 
#  Romberg : 21.333333333333332
# x, trapèzes : 22.0 
#  Simpson : 21.333333333333332

La fonction integrate.romb() n'admet pas de valeurs de x, elle ne fonctionne qu'avec un pas dx constant.

Notez que pour toutes les méthodes ci-dessus, si y est une matrice carrée, on peut indiquer l'axe d'intégration avec le paramètre axis.

Produit de convolution

Plusieurs fonctions permettent de faire un produit de convolution, mais il s'agit ici de la notion de filtre linéaire en traitement du signal (convolution discrète) : si l'on a deux vecteurs (a)_{0, M – 1} et (b)_{0, N – 1}, le produit de convolution vaut :

(a*b)_{n}=\sum _{m=-\infty }^{+\infty }a_{n-m}\cdot b_{m}

Par rapport au produit de convolution :

(f*g)(x)=\int _{t=-\infty }^{+\infty }f(x-t)\cdot g(t)\cdot \mathrm {d} t

la convolution discrète est une approximation à la largeur dt près. Pour avoir une approximation du produit de convolution, il faut donc corriger de ce facteur.

Pour l'échelle en x : le nombre de points calculé vaut N + M – 1, ce qui correspond aux différentes superpositions des vecteurs. Par exemple, pour un vecteur a de dimension 2 et un vecteur b de dimension 3 :

(a*b)_{0}=a_{0}b_{0}

(pour m > 0, le vecteur a_{n – m} n'est pas défini) ;

(a*b)_{1}=a_{1}b_{0}+a_{0}b_{1}

;

(a*b)_{2}=a_{1}b_{1}+a_{0}b_{2}

;

(a*b)_{3}=a_{1}b_{2}

.

Nous représentons les valeurs des vecteurs par des barres verticales pour représenter les « fenêtres glissantes » :

    |-|     |-|   |-|     |-|
|-|-|     |-|-|   |-|-|     |-|-|

nous avons bien 4 (2 + 3 – 1) états de superpositions.

Si le vecteur a correspond à la fonction ƒ définie sur [x_ƒ1 ; x_ƒ2] et le vecteur b à la fonction g définie sur [x_g1 ; x_g2], alors le point le plus à gauche (premier état de superposition) correspond à x = (x_ƒ1 – x_g2) ; le point le plus à droite (dernier état de superposition) correspond à x = (x_ƒ2 – x_g1). Il faut donc aussi corriger l'échelle en x.

Produit de convolution d'une fonction porte et d'une fonction demi-cercle.

La fonction la plus immédiate pour faire un produit de convolution est numpy.convolve(), avec en entrée deux vecteurs. Par exemple :

import numpy as np
import matplotlib.pyplot as plt

n = 100 # nombre de points
xmin = -2 # limites du calcul
xmax = 2
x = np.linspace(xmin, xmax, n)
dt = (xmax - xmin)/n # pas

def cercle(x):
    """Fonction cercle (équation cartésienne)"""
    bool = np.logical_and(x >= -1, x <= 1)
    y = np.zeros_like(x)
    y[bool] = np.sqrt(1 - x[bool]*x[bool])
    return y

def porte(x):
    """Fonction porte"""
    bool = np.logical_and(x >= -1, x <= 1)
    y = np.zeros_like(x)
    y[bool] = 1
    return y

y1 = cercle(x)
y2 = porte(x)

# Tracé des fonctions
plt.figure()
plt.plot(x, y1)
plt.plot(x, y2)
plt.gca().set_aspect("equal", adjustable="box")

# Calcul du produit de convolution
X = np.linspace(xmin - xmax, xmax - xmin, 2*n - 1)
Y = dt*np.convolve(y1, y2)

# Tracé du produit de convolution
plt.figure()
plt.plot(X, Y)

print(Y.shape, Y.max())
# (39,) 1.5159092533763792

Nous avons bien un maximum valant à peu près π/2 (aire du demi-disque de rayon 1).

Nous pouvons aussi détourner la fonction de multiplication des polynômes. Considérons les polynômes (∑a_i·xⁱ) et (∑b_i·xⁱ). Le produit de ces deux polynômes donne le polynôme :

(a_{0}+a_{1}\cdot x)\times (b_{0}+b_{1}\cdot x+b_{2}\cdot x^{2})=a_{0}b_{0}+(a_{0}b_{1}+a_{1}b_{0})x+(a_{0}b_{2}+a_{1}b_{1})x^{2}+a_{1}b_{2}x^{3}

Les coefficients de ce polynôme sont les mêmes que la convolution discrète. On peut donc utiliser la fonction de produit de polynômes, numpy.polynum() (ou numpy.polynomial.polynomial.polymul()) pour faire une convolution discrète.

L'extension scipy fournit plusieurs autres fonctions, dont certaines utilisent la transformée de Fourier rapide (FFT : fast Fourier-transform) : en effet, le produit de convolution correspond à un produit dans l'espace de Fourier. On a donc :

scipy.signal.fftconvolve() : utilise la transformée de Fourier ;
scipy.signal.convolve()() : choisit entre la méthode directe (somme) et la méthode de la transformée de Fourier, selon ce qui est le plus rapide ;
scipy.signal.oaconvolve() : utilise la méthode d'empiètement additif (overlap-add method).

Dans tous les cas, les vecteurs doivent avoir la même taille.

Équations différentielles

Équations différentielles ordinaires d'ordre 1

Une équation différentielle ordinaire (ÉDO) d'ordre 1 est une équation de la forme :

y'(t)=f(y(t))

que l'on écrit souvent sous la forme :

{\frac {\mathrm {d} y}{\mathrm {d} t}}=f(y)

Étant données des conditions initiales y(t₀) = y₀, on cherche la valeur de y pour un paramètre t_f.

Le module scipy contient un module integrate qui propose la fonction solve_ivp() ; c'est une procédure itérative. La syntaxe est de la forme scipy.integrate.solve_ivp(f, [t0, tf], y0). Le résultat est un dictionnaire contenant notamment les attributs :

.t : vecteur contenant les valeurs de t utilisées au cours de la procédure itérative ; la première valeur .t[0] est t₀, la dernière valeur .t[-1:] est la valeur t_f ;
.y : vecteur contenant les valeurs de la fonction ; la dernière valeur .y[-1:] est la valeur cherchée y(t_f), les autres valeurs sont celles évaluées au cours des itérations ;
.success : booléen valant True si le calcul a convergé, False s'il a échoué.

Par exemple, l'ÉDO d'ordre 1 :

y'(t)=t

a pour solution

y(t)={\frac {1}{2}}(t-t_{0})^{2}+y_{0}.

Nous considérons ici y(t₀) = 0 et t₀ = 0 donc y(t ) = ½t². Le code pour obtenir la solution est le suivant :

from scipy import integrate

def f(t, y):
    return t

y0 = [0]
t0 = 0
tf = 2
solu = integrate.solve_ivp(f, [t0, tf], y0, vectorized=True)

print("tf =", solu.t[-1:][0], "; y =", solu.y[0, -1:][0],
     "\n", "0,5⋅tf² =", 0.5*tf*tf)
# tf = 2.0 ; y = 2.0
# 0,5⋅tf² = 2.0

plt.plot(solu.t, solu.y.flat)
plt.plot(solu.t, solu.y.flat, "+")

Si la fonction est écrite de manière vectorisée, c'est-à-dire qu'elle peut s'appliquer à une matrice, on a alors intérêt à ajouter le paramètre vectorized=True dans l'appel de la fonction integrate.solve_ivp() : le calcul de la jacobienne est plus rapide.

Autres paramètres

Pour calculer la valeur à l'instant t_f, le solveur évalue successivement la valeur de la fonction à différents instants. On peut donc indiquer des valeurs intermédiaires avec le paramètre t_eval = [t0, t1, …, tn] ; le solveur renvoie alors les valeurs de la fonction à ces instants ce qui permet par exemple de tracer la courbe y(t ).

La fonction integrate.solve_ivp() peut utiliser plusieurs algorithmes. L'algorithme par défaut est la méthode de Runge-Kutta d'ordre 4 mais en utilisant l'ordre 5 à chaque étape ; cette méthode est appelée RK45. On peut choisir d'autre méthodes avec le paramètre method : RK23 pour une méthode de Runge-Kutta d'ordre 2, Radau, BDF (backward differentiation formula) et LSODA. De manière générale, on distingue :

les solveurs explicites, adaptés aux équations peu raides : RK45, RK23 ;
les solveurs implicites, adaptés aux équations raides : Radau et BDF ;
la méthode LSODA évalue la raideur et choisit soit algorithme Adams pour les problèmes non-raides, soit le BDF pour les problèmes raides.

Pour plus de détails voir : w:fr:Équation différentielle raide et w:fr:Méthodes de Runge-Kutta.

Enfin, le solveur integrate.solve_ivp() ne propose pas de manipuler une fonction ayant des paramètres supplémentaires, de type f(t, y, p1, p2, …). Pour cela, il faut encapsuler la fonction dans une autre fonction qui permet de fixer les paramètres :

def enveloppe(t, y):
    return f(t, y, p1, p2, …)

Il est possible pour cela d'utiliser la méthode lambda :

integrate.solve_ivp(lambda t, y: f(t, y, p1, p2), [t0, tf], y0)

Équations différentielles ordinaires d'ordres plus élevés

Une ÉDO d'ordre supérieur peut se réduire à une ÉDO d'ordre 1.

Pour plus de détails voir : w:fr:Équation différentielle ordinaire#Réduction à 1 de l'ordre d'une équation.

Nous pouvons donc utiliser la même fonction grâce à un changement de variables. Considérons par exemple l'ÉDO linéaire à coefficients constants d'ordre 2 suivante :

y ’’ + ay ’ + by + c = 0.

En posant

x₁ = y

x₂ = y ’

l'équation s'écrit

x₂’ + ax₁’ + bx₁ + c = 0.

En posant le vecteur

\mathrm {Y} ={\begin{pmatrix}x_{1}\\x_{2}\end{pmatrix}}

on a

\mathrm {Y} '={\begin{pmatrix}x_{1}'\\x_{2}'\end{pmatrix}}={\begin{pmatrix}x_{2}\\-ax_{1}'-bx_{1}-c\end{pmatrix}}.

En physique, on s'intéresse souvent à des fonctions du temps et l'on parle d'oscillations harmoniques ; on écrit souvent cette équation sous la forme :

{\ddot {s}}+2\lambda {\dot {s}}+\omega _{0}^{2}s+c=0

où λ est un terme d'amortissement et ω₀ la pulsation périodique ou pseudo-périodique. On a donc :

{\dot {\mathrm {Y} }}={\begin{pmatrix}x_{2}\\-2\lambda {\dot {x}}_{1}-\omega _{0}^{2}x_{1}-c\end{pmatrix}}.

On définit également le facteur de qualité Q :

\mathrm {Q} ={\frac {\omega _{0}}{2\lambda }}.

Pour plus de détails voir : w:fr:Oscillateur harmonique, w:fr:Équation différentielle linéaire d'ordre deux et v:fr:Équation différentielle/Équation différentielle linéaire du deuxième ordre à coefficients constants.

Prenons par exemple le cas d'une masse suspendue à un ressort. Le principe fondamental de la dynamique (PFD) s'écrit :

m{\ddot {x}}+kx+p=0

où :

m est la masse en kilogrammes (kg) ;
k est la raideur du ressort en newton par mètre (N/m) ;
p est le poids de l'objet, p = mg, l'accélération de la gravité valant g ≃ 9,81 m/s² ;
$x$ est l'altitude de l'objet (le ressort est détendu lorsque x = 0), ${\dot {x}}$ est la vitesse (en mètre par seconde, m/s) et ${\ddot {x}}$ l'accélération (en mètre par seconde au carré, m/s²) ;

soit

{\ddot {x}}+{\frac {k}{m}}x+{\frac {p}{m}}=0\Longleftrightarrow {\ddot {x}}+\omega _{0}^{2}x+g=0

en posant la pulsation, exprimée en s^–1 :

\omega _{0}={\sqrt {\frac {k}{m}}}.

Nous avons donc ici :

$x_{1}=x$ , la position ;
$x_{2}={\dot {x}}$ , la vitesse ;
${\dot {\mathrm {Y} }}={\begin{pmatrix}x_{2}\\-\omega _{0}^{2}{\dot {x}}_{1}-g\end{pmatrix}}.$

Nous prenons arbitrairement une pulsation ω = 1 s^–1, et comme conditions initiales x = 0 et x ’ = 0 à t = 0.

from scipy import integrate

k = 1 # raideur du ressort en N/m
m = 1 # masse du corps en kg
lmbd = 0
omega02 = np.sqrt(k/m)
g = 9.81

def f(t, Y):
    return [Y[1], -2*lmbd*Y[1] - omega02*Y[0] - g]

Y0 = [0, 0]
t0 = 0
n = 10
t = np.linspace(0, 2*np.pi, n)

solu = integrate.solve_ivp(lambda t, y: f(t, y, lmbd, omega02),
                           [t0, t[i]], Y0)
pos = solu.y[0, :]
vit = solu.y[1, :]

plt.plot(t, pos, label="position")
plt.plot(t, vit, label="vitesse")
plt.legend()

On retrouve bien des courbes sinusoïdales.

Considérons maintenant un circuit électrique RLC en série soumis à un échelon de tension E à t = 0. La loi des mailles donne

\mathrm {E} =\mathrm {L} {\frac {\mathrm {d} i}{\mathrm {d} t}}+\mathrm {R} i+{\frac {q}{\mathrm {C} }}

soit

{\ddot {q}}+{\frac {\mathrm {R} }{\mathrm {L} }}{\dot {q}}+{\frac {q}{\mathrm {L} \mathrm {C} }}={\frac {\mathrm {E} }{\mathrm {L} }}

donc

\lambda ={\frac {\mathrm {R} }{2\mathrm {L} }}\ ;\ \omega _{0}^{2}={\frac {1}{\mathrm {L} \mathrm {C} }}.

from scipy import integrate

L = 50 # inductance en H
R = 200 # résistance en Ω
C = 1e-5 # capacité en F
E = 1 # échelon de tension en V

foo = 1/(L*C)
lmbd = 0.5*R*C*foo
omega02 = foo

def equation(t, y, lmbd, omega02):
    return [y[1], -2*lmbd*y[1] - omega02*y[0] + E/L]

t0 = 0

T = 2*np.pi/np.sqrt(omega02 - lmbd*lmbd)
tmax = 5/lmbd
t = np.arange(0, tmax, T/8)
n = len(t)

y0 = [0, 0]

resultat = integrate.solve_ivp(lambda t, y: equation(t, y, lmbd, omega02),
                                   [t0, t[-1:]], y0, t_eval=t, method="Radau")
q = resultat.y[0, :]
i = resultat.y[1, :]

fig, axes1 = plt.subplots(1, 1)
axes2 = axes1.twinx()
axes1.plot(t, q, "orange", label="charge")
axes1.legend(loc=1)
axes1.set_xlabel("$t$ (s)")
axes1.set_ylabel("$q$ (C)")
axes2.plot(t, i, label="intensité")
axes2.legend(loc=4)
axes2.set_ylabel("$i$ (A)")

On retrouve bien une courbe de la forme $q(t)=\mathrm {A} e^{-\lambda t}\cos \left(\omega _{0}{\sqrt {1-{\frac {\lambda ^{2}}{\omega _{0}^{2}}}}}+\varphi \right).$

SUNDIALS

Le LLNL (Lawrence Livemore National Laboratory) a développé une bibliothèque pour le calcul différentiel et la résolution de systèmes algébriques, appelé SUNDIALS (Suite of nonlinear and differential/algebraic equation solvers)^[1]. Cette bibliothèque est actuellement en cours de transposition dans Python^[2].

Notes et références

↑ « SUNDIALS: SUite of Nonlinear and DIfferential/ALgebraic Equation Solvers », sur LLNL (consulté le 21 mars 2023).
↑ « python-sundial », sur PyPi (consulté le 21 mars 2023) ; ne pas confondre avec le module sundial (au singulier) qui sert à afficher des barres de progression (le mot sundial signifiant « cadran solaire ».

Algèbre linéaire < ↑ > Résolution d'équations

[1] « SUNDIALS: SUite of Nonlinear and DIfferential/ALgebraic Equation Solvers », sur LLNL (consulté le 21 mars 2023).

[2] « python-sundial », sur PyPi (consulté le 21 mars 2023) ; ne pas confondre avec le module sundial (au singulier) qui sert à afficher des barres de progression (le mot sundial signifiant « cadran solaire ».

[1]

[2]