Python pour le calcul scientifique/Résolution d'équations

Rappelons que dorénavant les programmes commencent tous par :

#!/usr/bin/python3

import numpy as np
import matplotlib.pyplot as plt

Système d'équations linéaires

Le sujet a été abordé dans le chapitre concernant l'algèbre linéaire. Modèle:Loupe

Équation polynomiale

La résolution d'une équation polynomiale consiste à trouver les racines de son polynôme. Par exemple, pour résoudre l'équation

$x^2+3x+5=0$

nous pouvons utiliser :

import numpy.polynomial.polynomial as nppol

X = nppol.polyroots([5, 3, 1]) # [-1.5-1.6583124j -1.5+1.6583124j]

Équation différentielle

Le sujet est abordé dans le chapitre sur le calcul différentiel. Modèle:Loupe

Équation quelconque

Une équation peut s'écrire de manière générale

y = ƒ(x)

ce qui revient à résoudre

ƒ(x) – y = 0.

Le module scipy.optimize propose la fonction root() pour cela.

Par exemple, nous cherchons ci-dessous à résoudre x² + 5⋅sin x = 2

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize as opt

def f(x):
    '''fonction présentant au moins un zéro'''
    return x*x + 5*np.sin(x)

def residu(x, y):
    return f(x) - y

x = np.linspace(-4, 4, 40)
y = f(x)

ycible = 2
x0 = 0 # point de départ

minimum = opt.root(residu, x0, args=(ycible))
xopt = minimum.x[0]
yopt = f(xopt)
print(f"La valeur y = {round(yopt, 2)} est atteinte pour x = {round(xopt, 2)}")
# La valeur y = 2.0 est atteinte pour x = 0.38

plt.plot(x, y)
plt.plot((x[0], x[-1]), (ycible, ycible), "--k", linewidth=0.5)
plt.plot(xopt, yopt, "+k")

Le graphique nous montre qu'il existe une autre solution (x = –2,35). Nous pouvons trouver la trouver avec un point de départ différent, par exemple x0 = -5.

Cette méthode marche aussi pour des systèmes d'équations, en recherchant le vecteur nul. Considérons par exemple une poutre de section rectangulaire L × l, par exemple une planche de largeur L et de hauteur l. Sa raideur dépend du matériau dont elle est faite (bois, acier, aluminium…) mais aussi de sa forme ; la contribution de la forme à la raideur est exprimée par le moment quadratique I, exprimé en mm⁴ (si L et l sont en mm). Mise à plat, le moment quadratique vaut

I₁ = L × l ³/12.

Mise à chant (c'est-à-dire posée sur la tranche), il vaut :

I₁ = L³ × l /12.

Nous cherchons les dimensions L et l telles que I₁ = Modèle:Unité et I₂ = Modèle:Unité. Pour cela, nous définissons le vecteur d'état X = ([L, l]) et la fonction moment quadratique I(X) = ([I1, I2]). Nous partons de dimensions a priori L = l = Modèle:Unité. Nous cherchons donc à résoudre le système

$\{\begin{array}{c}\mathrm{L}\times l^3=200000\\ {\mathrm{L}}^3\times l=700000\end{array}$

par l'équation vectorielle

$\mathrm{I}\left(\begin{array}{c}\mathrm{L}\\ l\end{array}\right)-\left(\begin{array}{c}200000\\ 700000\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right).$

import numpy as np
from scipy import optimize as opt

def I(X):
    L = X[0]
    LLL = L*L*L # L³
    l = X[1]
    lll = l*l*l # l³
    return (1/12)*np.array([L*lll, LLL*l]) # ([I1, I2]) en mm⁴

def residu(X, Y):
    return I(X)-Y

ycible = np.array([200000, 700000])

x0 = np.array([50, 50])

minimum = opt.root(residu, x0, args=(ycible))

print("Poutre\n")
print(minimum)
print("\n")
Xopt = minimum.x
Lopt = Xopt[0]
lopt = Xopt[1]
print(f"Section : ({Lopt.:3f}, {lopt:.3f}) mm") # (62.962, 33.655) mm
print(f"Moments quadratiques : ({I(Xopt)}) mm⁴") # ([199999.99999985 700000.00000009]) mm⁴

Il nous faut donc une planche de dimensions approximatives 63 × 34 mm².

Notes et références

---

[[../Calcul différentiel et intégral|Calcul différentiel et intégral]] < [[../|↑]] > [[../Calcul symbolique|Calcul symbolique]]