Calcul scientifique sous Python

Introduction #

Dans ce chapitre, nous nous intéressons à des méthodes numériques pour la résolution d'équations sur les réels, d'équations différentielles et le calcul approché d'intégrales. Nous rappelons tout d'abord quelques fonctions mathématiques qui sont faciles à programmer, puis nous expliquons comment utiliser les fonctions prédéfinies de Python et ses bibliothèques dédiées pour ce type de calcul.

Les fonctions présentées dans ce chapitre nécessitent l'importation des modules suivants :

• import numpy as np
• import scipy.optimize as resol
• import scipy.integrate as integr
• import matplotlib.pyplot as plt

Avec :

• NumPy : Le module numpy est utilisé pour manipuler des vecteurs et des matrices de nombres.
• SciPy : Le module scipy regroupe un certain nombre de sous modules qui sont très riches en fonctions prédéfinies pour assurer plusieurs routines. Scipy est basé sur numpy, mais il en étend considérablement les possibilités de ce dernier et il permet de faire des calculs en statistiques, optimisation, intégration numérique, traitement du signal, traitement d'image, algorithmes génétiques, etc.
• MatPlotLib : Matplotlib est une bibliothèque en Python très utilisée pour tracer des graphiques en deux et trois dimensions.

Le module NumPy #

Le module numpy permet d'effectuer des calculs sur des vecteurs ou des matrices, élément par élément, via un nouveau type d'objet appelé array. Ce module contient plusieurs fonctions pour faire par exemple de l'algèbre linéaire et la génération de nombre aléatoire.

Voici des exemples d'utilisation de plusieurs fonctions du module NumPy :

• array(L) : crée un tableau à partir d'une liste L.

>>> import numpy as np
>>> L=[1,2,3]
>>> T=np.array(L) # créer un tableau à partir d'une liste
>>> type(T)
<class 'numpy.ndarray'>
>>> T
array([1, 2, 3])

Il est également possible de construire des objets array à deux dimensions, il suffit de passer en argument une liste de listes à la fonction array().

>>> np.array([[1,2],[3,4]]) # Créer une matrice à partir d'une liste de listes
array([[1, 2],
       [3, 4]])

• arange(a, b, k) : est équivalente à range() et permet de construire un array à une dimension de manière simple.

>>> np.arange(10)
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> np.arange(2,10)
array([2, 3, 4, 5, 6, 7, 8, 9])
>>> np.arange(2,10,2)
array([2, 4, 6, 8])
>>> np.arange(10,2,-2)
array([10,  8,  6,  4])

Un des avantages de la fonction arange() est qu'elle permet de générer des objets array qui contiennent des entiers ou des réels selon l'argument qu'on lui passe.

>>> np.arange(2,3,0.1)
array([ 2. ,  2.1,  2.2,  2.3,  2.4,  2.5,  2.6,  2.7,  2.8,  2.9])

Pour récupérer un ou plusieurs élément(s) d'un objet array, vous pouvez utiliser l'indiçage ou les tranchages (slicing), de la même manière que pour les listes.

>>> T=np.arange(2,10)
>>> T
array([2, 3, 4, 5, 6, 7, 8, 9])
>>> T[0] # élément d'indice 0
2
>>> T[1:4] # les éléments d'indices 0,1,2 et 3
array([3, 4, 5])
>>> T[5:] # du 5ème au dernier élément
array([7, 8, 9])
>>> T[:5] # du début au 4ème élément
array([2, 3, 4, 5, 6])
>>> T[::2] # tous les éléments d'indice pair
array([2, 4, 6, 8])

Dans le cas d'un objet array à deux dimensions, vous pouvez récupérer une ligne, une colonne ou bien un seul élément.

>>> T=np.array([[1,2,3],[4,5,6]])
>>> T
array([[1, 2, 3],
       [4, 5, 6]])
>>> T[0,0] # élément d'indice (0,0)
1
>>> T[0,:] # une ligne entière, ligne numéro 0
array([1, 2, 3])
>>> T[0] # une ligne entière, ligne numéro 0
array([1, 2, 3])
>>> T[:,2] # une colonne entière, colonne numéro 2
array([3, 6])
>>> T[[0,2]] # les lignes numéros 0 et 2
array([[1, 2, 3],
       [7, 8, 9]])

• linspace(a, b, n) : crée un vecteur de n valeurs régulièrement espacées entre a et b (inclus).

>>> T=np.linspace(0,10,5)
>>> T
array([  0. ,   2.5,   5. ,   7.5,  10. ])

• zeros(n) : crée un tableau de taille n rempli de zéros.
• zeros((n, m)) : crée une matrice de taille (n, m) rempli de zéros.

>>> z=np.zeros((2,2)) # matrice nulle de taille 2*2
>>> z
array([[ 0.,  0.],
       [ 0.,  0.]])
>>> b=np.zeros((2,2),dtype=np.bool) # changer le type des éléments de la matrice
>>> b
array([[False, False],
       [False, False]], dtype=bool)

• ones(n) : crée un tableau de taille n rempli de uns.
• ones((n, m)) : crée une matrice de taille (n, m) rempli de uns.

>>> o=np.ones((2,3),int) # créer une matrice d'entiers contenant que des uns
>>> o
array([[1, 1, 1],
       [1, 1, 1]])

Par défaut, les fonctions zeros() et ones() génèrent des réels, mais vous pouvez demander des entiers en passant l'option int en second argument.

• eye(n) : crée une matrice identité de taille (n, n).

>>> I=np.eye(3) # matrice identité 3*3
>>> I
array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])

• diag(L) : crée une matrice diagonale.

>>> d=np.diag([1,2,3])
>>> d
array([[1, 0, 0],
       [0, 2, 0],
       [0, 0, 3]])

• concatenate((a₁, a₂, ⋯, a_n), axis=1 ou 0) : permet d'accoler deux ou plusieurs tableaux (horizontalement avec axis=1, verticalement avec axis=0). En cas de concaténation horizontale par exemple, toutes les colonnes doivent avoir la même hauteur. Les différents tableaux à concaténer doivent être donnés sous la forme d'un tuple.

>>> A=np.random.randint(1,6,(2,3))
array([[4, 3, 2],
       [4, 1, 2]])
>>> B=np.random.randint(7,11,(2,4))
array([[10, 7, 10, 9],
       [ 9, 8,  7, 9]])
>>> np.concatenate((A,B),axis=1)
array([[ 4, 3, 2, 10, 7, 10, 9],
       [ 4, 1, 2,  9, 8,  7, 9]])

• T.shape : pour obtenir la taille du tableau T.
• T.dtype : indique le type des éléments de T.
• T.size : indique le nombre d'éléments du tableau T.
• T.ndim : le nombre d'indices nécessaires au parcours du tableau, c'est-à-dire le nombre d'éléments du tuple indiquant son format.

>>> T=np.array([[1,2,3],[4,5,6]])
>>> T.shape # un tuple (nombre de ligne, nombre de colonne)
(2, 3)
>>> T.dtype
dtype('int32')
>>> T.size
6
>>> T.ndim
2

• reshape(M, (dimensions)) et resize(M, (dimensions)) : les fonctions reshape() et resize() permettent de redimensionner à volonté les dimensions d'un array. Il faut leur passer en argument l'objet array à redimensionner ainsi qu'un tuple indiquant la nouvelle dimension.

>>> T=np.arange(8)
>>> a=np.reshape(T,(2,4)) # redimensionner le vecteur en matrice 2*4
>>> a
array([[0, 1, 2, 3],
       [4, 5, 6, 7]])
>>> b=np.reshape(T,(4,2))
>>> b
array([[0, 1],
       [2, 3],
       [4, 5],
       [6, 7]])
>>> c=np.reshape(T,(4,3))
ValueError: total size of new array must be unchanged

La fonction reshape() attend un tuple dont la dimension est compatible avec le nombre d'éléments contenus dans l'array de départ, alors que resize() remplira le nouvel objet array généré même si les longueurs ne coincident pas.

>>> c=np.resize(T,(4,3))
>>> c
array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 0],
       [1, 2, 3]])

• dot(M, N) : pour effectuer un produit matriciel de deux matrices.

>>> A=np.array([[1,2],[3,4]])
>>> A
array([[1, 2],
       [3, 4]])
>>> B=np.array([[1,1,1],[2,2,2]])
>>> np.dot(A,B) # produit matricielle de A par B
array([[ 5,  5,  5],
       [11, 11, 11]])

• vdot(v1, v2) : pour effectuer un produit scalaire de deux vecteurs.

>>> A=np.array([1,2])
>>> B=np.array([3,4])
>>> np.vdot(A,B) # le produit scalaire de A par B
11

• transpose(M) : pour transposer une matrice.

>>> T=np.array([[1,2,3],[4,5,6]])
>>> T
array([[1, 2, 3],
       [4, 5, 6]])
>>> np.transpose(T)
array([[1, 4],
       [2, 5],
       [3, 6]])

• rank(M) : rang d'une matrice.

>>> M=np.array([[1,2,3],[1,2,4]])
>>> np.rank(M)
2

• mean(T) : valeur moyenne d'un tableau.
• min(T) : Le minimum d'un tableau.
• max(T) : Le maximum d'un tableau.
• sum(T) : La somme des valeurs d'un tableau.

>>> M=np.array([[1,2,3],[1,2,4]])
>>> np.mean(M)
2.1666666666666665
>>> np.min(M)
1
>>> np.max(M)
4
>>> np.sum(M)
13

• linalg.inv(M) : inversion d'une matrice.
• linalg.det(M) : déterminant d'une matrice.

>>> M=np.array([[1,2],[3,4]])
>>> np.linalg.det(M)
-2.0000000000000004
>>> np.linalg.inv(M)
array([[-2. ,  1. ],
       [ 1.5, -0.5]])

• linalg.eig(M) : renvoie un tuple dont le premier élément correspond aux valeurs propres et le second élément aux vecteurs propres.

>>> M=np.array([[1,2],[3,4]])
>>> np.linalg.eig(M)
(array([-0.37228132,  5.37228132]), array([[-0.82456484, -0.41597356],
       [ 0.56576746, -0.90937671]]))
>>> np.linalg.eig(M)[0]
array([-0.37228132,  5.37228132])
>>> np.linalg.eig(M)[1]
array([[-0.82456484, -0.41597356],
       [ 0.56576746, -0.90937671]])

• linalg.solve(A, b) : résolution du système linéaire AX = b.

>>> A,b=np.array([[1,2],[3,4]]), np.array([[1,1]])
>>> x=np.linalg.solve(A,np.transpose(b))
>>> np.dot(A,x)
array([[ 1.],
       [ 1.]])

Les opérations sur des tableaux

Les opérations +, *, /, ==, etc. s'appliquent aux tableaux numpy, mais terme à terme.

>>> A=np.array([2,2,2])
>>> B=np.array([3,3,3])
>>> A+B
array([5, 5, 5])
>>> A-B
array([-1, -1, -1])
>>> A**B
array([8, 8, 8], dtype=int32)
>>> A/B
array([ 0.66666667,  0.66666667,  0.66666667])
>>> A//B
array([0, 0, 0], dtype=int32)
>>> A%B
array([2, 2, 2], dtype=int32)
>>> A*B
array([6, 6, 6])
>>> A<B
array([ True,  True,  True], dtype=bool)
>>> A!=B
array([ True,  True,  True], dtype=bool)
>>> A==B
array([False, False, False], dtype=bool)

Les fonctions universelles appliquées sur des tableaux

Les fonctions définies dans le module numpy sont universelles, ou vectorialisées, c'est-à-dire qu'elles acceptent comme argument un tableau numpy et renvoient le tableau numpy des images.

>>> np.exp(A)
array([  2.71828183,   7.3890561 ,  20.08553692])
>>> np.log(A)
array([ 0.        ,  0.69314718,  1.09861229])
>>> f=lambda x:x**2
>>> f(2)
4
>>> f(A)
array([1, 4, 9])

Les tableaux pseudo-aléatoires

Les fonctions qui forment des tableaux pseudo-aléatoires sont présentes dans le sous-module random du module numpy.

• random.rand(n) : renvoie un tableau de n valeurs pseudo-aléatoires et uniformément distribuées dans l'intervalle [0, 1].

>>> import numpy as np
>>> x=np.random.rand(6)
>>> x
array([ 0.9595917, 0.92122968, 0.3544391,
        0.10543147, 0.81636296,
        0.44601582])

• random.randn(n) : suit la même syntaxe que random.rand, mais elle renvoie un échantillon de n valeurs pseudo-aléatoires au sens de la loi normale centrée réduite (c'est-à-dire d'espérance 0 et d'écart-type 1).

import matplotlib.pyplot as plt
import numpy as np
# 1000 tirages de loi normale de moyenne 0 et d'écart type 1
x = np.random.randn(1000) # Génération de n valeurs de loi normale centrée réduite
n, bins, patches = plt.hist(x, 50, normed=1, facecolor='b', alpha=0.5) # La valeur 50 donné en deuxième paramètre de la méthode hist indique le nombre de barres à afficher.
plt.xlabel('x')
plt.ylabel('Probabilité de x N(0,1)(x)')
plt.axis([-5, 5, 0, 0.5])
plt.grid(True)
plt.show()

• random.randint(a, b) : renvoie une valeur entière pseudo-aléatoire au sens de la distribution uniforme dans un intervalle semi-ouvert [a, b[.
• random.randint(a, b, (n, m)) : renvoie un tableau de taille n × m valeurs entières pseudo-aléatoires au sens de la distribution uniforme dans un intervalle semi-ouvert [a, b[. Si on omet la valeur a, les entiers pseudo-aléatoires sont choisis dans [0, b[.

>>> np.random.randint(5)
1
>>> np.random.randint(1,7) # un lancer de dé
5
>>> np.random.randint(1,7,3) # trois lancers de dé
array([2, 5, 4])
>>> np.random.randint(1,7,(3,2)) # 6 lancers de dé sous forme d'une matrice 3x2
array([[1, 1],
       [4, 5],
       [6, 1]])

import matplotlib.pyplot as plt
import numpy as np
# 1000 tirages entre 1 et 6
x = np.random.randint(1,7,1000)
n, bins, patches = plt.hist(x, 6, normed=1, facecolor='pink', alpha=0.8)
plt.xlabel('dé')
plt.ylabel('Probabilité')
plt.axis([1, 6, 0, 0.3])
plt.grid(True)
plt.show()

Le module MatPlotLib #

Pour le simple tracé de courbes nous n'utiliserons que le sous-module pyplot, importé, avec alias, à l'aide de l'instruction : import matplotlib.pyplot as plt. Pour plus de documentation : http://www.matplotlib.org.

Les fonctions essentielles de pyplot sont :

• plot() pour le tracé de points, de courbes
• Les fonctions xlabel() et ylabel() sont utiles pour donner un nom aux axes.
• title() permet de définir le titre du graphique.
• grid() affiche une grille en filigrane.
• show() pour afficher le graphique créé.
• savefig(fich) : permet de sauvegarder la figure en cours dans le fichier fich.

Utiliser plot() avec :

• en 1^er argument la liste des abscisses,
• en 2^e argument la liste des ordonnées,
• en 3^e argument (optionnel) le motif des points :
  • '.' pour un petit point,
  • 'o' pour un gros point,
  • '+' pour une croix,
  • '*' pour une étoile,
  • '-' points reliés par des segments,
  • '--' points reliés par des segments en pointillés,
  • '-o' gros points reliés par des segments (on peut combiner les options),
  • 'b', 'r', 'g', 'y' pour de la couleur (bleu, rouge, vert, jaune, etc.).

Recherche approchée de f(x) = 0 #

Dans les problèmes d'ingénierie, il faut souvent rechercher les solutions d'une équation de la forme f(x) = 0, où f est une fonction à valeurs réelles. Le but est donc de trouver une solution approchée tout en minimisant le temps d'exécution et en obtenant une bonne approximation. Deux méthodes classiques de recherche de zéro sont présentées dans ce cours, à savoir :

• La méthode de dichotomie
• La méthode de Newton

On considérera par la suite une fonction f continue sur un intervalle [a, b] à valeurs réelles, telle que f(a) · f(b) < 0. D'après le théorème des valeurs intermédiaires, la fonction f s'annule au moins une fois sur l'intervalle [a, b].

La méthode de dichotomie

Écrire une fonction Python dichotomie(f, a, b, e) qui recherche et renvoie un zéro de f dans l'intervalle [a, b] et selon la précision e.

Code Python pour la méthode par dichotomie :

def dichotomie(f,a,b,e):
    while b-a > e:
        m = (a+b)/2
        if f(a)*f(m) < 0:
            b = m # Dichotomie à gauche
        else:
            a = m # Dichotomie à droite
    return m

def f(x):
    return x**2-2

print(dichotomie(f,0,2,10**(-6)))
# résultat 1.4142141342163086

from math import *
print(sqrt(2))
# résultat 1.4142135623730951

La méthode de Newton

Écrire une fonction Python newton(f, g, u0, e) qui recherche et renvoie un zéro de f selon la précision e, avec g la fonction dérivée de f et u0 la valeur initiale.

Code Python pour la méthode de Newton :

def newton(f, g, u0, e):
    u = u0
    # la fonction g est la dérivée de f
    v = u0 - f(u0)/g(u0)
    while abs(v-u) > e:
        u = v
        v = v - f(v)/g(v)
    return v

def f(x):
    return x**2-2

def g(x):
    return 2*x

print(newton(f,g,2,10**(-6)))
# résultat 1.4142135623730951

from math import *
print(sqrt(2))
# résultat 1.4142135623730951

Concernant les méthodes de résolution approchée d'équations, tout se trouve dans le module scipy.optimize, qui est lui-même un sous-module de scipy. Il contient entre autres les fonctions suivantes :

import scipy.optimize as resol

def f(x):
    return x**2-2

a=1
b=2
x=resol.bisect(f,a,b)
print("sol=",x)
x=resol.newton(f,a)
print("sol=",x)
x=resol.fsolve(f,a)
print("sol=",x)
x=resol.root(f,a)
print("sol=",sol.x)

Résoudre un système d'équations #

Résoudre graphiquement et numériquement le système d'équations suivant :

{ x − 2y = 3 ;    4x + 5y = 6 }

import matplotlib.pyplot as plt
import numpy as np

def d1(x):
    return (3-x)/2

def d2(x):
    return (6-4*x)/5

x=np.linspace(-2,2,100)
y1=d1(x)
y2=d2(x)
plt.plot(x,y1,'r')
plt.plot(x,y2,'b')
plt.grid("on")
plt.show()

Résoudre graphiquement et numériquement le même système d'équations.

La déclaration de la fonction f de la manière suivante :

def f(x):
    return [x[0]+2*x[1]-3, 4*x[0]+5*x[1]-6]

En considérant comme point de départ le vecteur (0, 0), on obtient :

print(resol.fsolve(f,(0,0)))
# Résultat: array([-1, 2])

Les fonctions vectorielles #

Une démarche strictement identique doit être suivie pour obtenir la solution à un système d'équations non linéaires, à condition de n'utiliser qu'un seul argument (list ou tuple) pour stocker les différentes composantes du vecteur x = (x₁, x₂, ⋯, x_n) considéré.

En guise d'exemple, considérons le système non-linéaire formé par les fonctions :

f₁(x₁, x₂) = x₁ + 3 · log₁₀(x₁) − x₂²

f₂(x₁, x₂) = 2 · x₁² − x₁ · x₂ − 5 · x₁ + 1

La déclaration de la fonction f est désormais classique :

def f(x):
    return [x[0]+3*log10(x[0])-x[1]**2, 2*x[0]**2-x[0]*x[1]-5*x[0]+1]

En considérant comme point de départ le vecteur, x₀ = (5.0, 5.0), on obtient :

print(resol.fsolve(f,(5.0,5.0)))
# Résultat: array([3.48744279, 2.26162863])

Les racines d'un polynôme #

Le module roots de la bibliothèque numpy détermine les racines dans ℂ d'un polynôme.

import numpy as np

P1=[1,0,-2]  # x**2-2
racines1=np.roots(P1)
print(racines1)
# résultat [ 1.41421356 -1.41421356]

P2=[1,0,1,1]  # x**3+x+1
racines2=np.roots(P2)
print(racines2)
# résultat [ 0.3411639+1.1615414j
#            0.3411639-1.1615414j
#           -0.6823278+0.j ]

Calcul approché d'une intégrale #

• Calculer une valeur approchée d'une intégrale n'est en général pas très difficile.
• Les méthodes les plus simples consistent à subdiviser régulièrement l'intervalle d'intégration pour approcher le calcul d'intégrale par un calcul d'aires de figures géométriques plus simples (par exemple dans la méthode des trapèzes).
• Subdivision régulière du segment [a; b] :
• Soit [a; b] un intervalle non vide et non réduit à un point (i.e. a < b).
• Une subdivision régulière de [a; b] en n+1 points est une suite finie de n+1 réels :
  • (a_k)_{0 ≤ k ≤ n} vérifiant :
    • a₀ = a, a_n = b
    • et pour tout k ∈ [0; n], a_k = a + k · (b − a)/n.
• Sous Python :

>>> import numpy as np
>>> x = np.linspace(a, b, n+1)

Méthode des rectangles

La méthode des rectangles revient à une approximation par une fonction en escalier, avec n « marches » de longueur (b−a)/n. La valeur approchée R de l'intégrale vaut alors :

R = ((b−a)/n) · Σ_i=0ⁿ⁻¹ f(a_i)

Exemple :

∫₀^π/2 sin(x) dx = 1

import numpy as np  # pour la fonction linspace()

def rectangle(f,a,b,n):
    x = np.linspace(a,b,n+1) # Subdivision régulière
    Y = f(x[:n]) # Séquence des f(ak)
    return (b-a)/n * sum(Y)


print(rectangle(np.sin,0,np.pi/2,10))
# résultat 0.919403170015
print(rectangle(np.sin,0,np.pi/2,1000))
# résultat 0.99921439622

Méthode des trapèzes

La valeur approchée R de l'intégrale vaut :

R = ((b−a)/n) · ((f(a) + f(b))/2 + Σ_i=1ⁿ⁻¹ f(a_i))

Exemple :

∫₀^π/2 sin(x) dx = 1

import numpy as np  # pour la fonction linspace()

def trapeze(f,a,b,n):
    x = np.linspace(a,b,n+1) # Subdivision régulière
    Y = f(x[1:n]) # Séquence des f(ak)
    return (b-a)/n * (sum(Y)+(f(a)+f(b))/2)


print(trapeze(np.sin,0,np.pi/2,10))
# résultat 0.997942986354
print(trapeze(np.sin,0,np.pi/2,1000))
# résultat 0.999999794383

scipy dispose de plusieurs méthodes d'intégration de fonctions dans son sous-module integrate :

>>> import scipy.integrate as integr

La méthode quad(f, a, b) pour intégrer f sur l'intervalle [a; b] :

>>> integr.quad(np.sin,0,np.pi/2)
0.9999999999999999, 1.1102230246251564e-14

Elle retourne un couple constitué de la valeur approchée de l'intégrale (1^er élément) et d'une estimation de l'erreur commise.

• trapz() applique la méthode des trapèzes à un échantillonnage :

>>> x = np.linspace(0,np.pi/2,1001)
>>> integr.trapz(np.sin(x),dx=(np.pi/2)/1000)
0.99999979438323383

• simps() applique la méthode de Simpson à un échantillonnage :

>>> x=np.linspace(0,np.pi/2,1001)
>>> integr.simps(np.sin(x),dx=(np.pi/2)/1000)
1.0000000000000338