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

Un livre de Wikilivres.
Sauter à la navigation Sauter à la recherche

Rappelons que dorénavant les programmes commencent tous par :

import numpy as np
import matplotlib.pyplot as plt

Dérivation numérique[modifier | modifier le wikicode]

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.

Enfin, la fonctopn 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.]

Équations différentielles[modifier | modifier le wikicode]

Équations différentielles ordinaires d'ordre 1[modifier | modifier le wikicode]

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

que l'on écrit souvent sous la forme :

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

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 t0, la dernière valeur .t[-1:] est la valeur tf ;
  • .y : vecteur contenant les valeurs de la fonction ; la dernière valeur .y[-1:] est la valeur cherchée y(tf), 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

Nous considérons ici y(t0) = 0 donc y(t ) = ½t 2. 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[modifier | modifier le wikicode]

Pour calculer la valeur à l'instant tf, 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[modifier | modifier le wikicode]

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

x1 = y
x2 = y ’

l'équation s'écrit

x2’ + ax1’ + bx1 + c = 0.

En posant le vecteur

on a

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 :

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

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


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 :

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/s2 ;
  • est l'altitude de l'objet (le ressort est détendu lorsque x = 0), est la vitesse (en mètre par seconde, m/s) et l'accélération (en mètre par seconde au carré, m/s2) ;

soit

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

Nous avons donc ici :

  • , la position ;
  • , la vitesse ;

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

soit

donc

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

Notes et références[modifier | modifier le wikicode]


Régression et optimisation < >