Calcul symbolique avec Python
Le calcul symbolique
Sur TangenteX, nous avons l'habitude de faire du calcul numérique. Il faut dire que c'est la première chose que l'on a fait faire à un ordinateur, qui s'appelle "computer" en anglais ! Par exemple, si nous cherchons la valeur de l'intégrale de f(x) = x entre 0 et 1, nous n'allons pas utiliser la primitive de f(x) mais plutôt un algorithme qui va calculer l'aire sous la courbe entre 0 et 1. Très satisfaisant en calcul numérique, mais parfois il serait utile de connaitre la valeur de cette primitive ! De même, en présence d'une expression très compliquée, nous ne demandrons pas à notre programme de la simplifions, mais nous prendrons la peine de la programmer ou bien de la simplifier "à la main".
Dans cette page, nous allons nous intéresser à une autre manière de calculer, celle que vous pratiquez le plus communément au lycée, en prépa ou à la fac : le calcul symbolique. Vous vous doutez peut-être que ce n'est pas la méthode naturelle de calcul d'un ordinateur ! Essayez d'imaginer un programme qui vous calcule la dérivée de sin(x). Tout d'abord, la question n'a pas vraiment de sens en calcul numérique, à défaut d'avoir mentionner en quel point je veux connaitre la dérivée. Puis, que peut bien signifier pour un ordinateur et ses programmes de calcul l'expression "sin(x)" ? Qu'est-ce que "x", quelle valeur lui attribuer ?
Il existe une famille de programmes, les CAS (Computer Algebra System) qui sont conçus pour effectuer de tels calculs, non plus avec des nombres uniquement, mais avec des symboles. Ils effectuent ce qu'on appelle du calcul formel. Vous en connaissez sans doute au moins un, Maple ou/et Mathematica, pratiqués en classe prépa et à la fac, très peu au lycée malheureusement. En logiciels libres, citons Xcas et Yacas que vous pourrez tester si le coeur vous en dit. Pour cette page, j'ai décidé d'utiliser le module sympy de Python, qui sans être parfait rend bien des services.
Pour réaliser cette page, je me suis très largement inspiré du tutorial du projet sympy.
Installer et initialiser sympy
L'environnement de travail
Comme d'habitude, j'utilise l'environnement Spyder2 de la suite Anaconda, en version 2.3.9. Vous pouvez aussi utiliser idle, mais il vous faudra charger une console iPython, car sympy s'utilise en mode interactif la plupart du temps. Spyder2 fonctionne sous Windows et OSX et j'utilise les deux versions indifféremment. Les manipulations suivantes sont faites sur une console iPython en Pyton 2.7.11.
Initialiser sympy
Le module sympy est en principe déjà installé avec Python. Pour vous en assurer, tapez l'instruction suivante sur la console iPython :
In [2]: from sympy import *
La chaine "In [2]:" est le prompt de la console iPython, le nombre entre crochets indiquant le rang de la commande passée. Comme vous l'avez deviné, cette instruction importe toutes les fonctions du module sympy. Si celui-ci est absent, la console vous affichera un message d'erreur et vous devrez charger ce module à la main. Mais c'est peu probable...
Nous allons travailler en interactif sur la console iPython et sympy va nous afficher des résultats. Pour qu'ils apparaissent de façon agréable, je vais initialiser la fonction d'affichage de la façon suivante :
In [3]: init_printing()
Essayez taper par exemple, l'instruction suivante (suivi de Return évidemment...)
In [4]: Integral(sqrt(1/x),x)
Vous verrez s'afficher votre intégrale en format LaTex du plus bel effet.
Vous pouvez également utiliser l'instruction
In [5]: init_session()
qui effectue l'intialisation de la librairie, de la visualisation mais aussi la déclaration des symboles les plus fréquemment utilisés :
from __future__ import division
from sympy import *
x, y, z, t = symbols('x y z t')
k, m, n = symbols('k m n', integer=True)
f, g, h = symbols('f g h', cls=Function)
init_printing()
Nous verrons ce qu cela signifie tout de suite. Pour ouvrir une session de calcul formel sous Python, je l'utilise systématiquement car elle nous simplifie bien la vie.
Cette initialisation étant effectuée, nous allons pouvoir passer aux choses sérieuses, et en premier déclarer des symboles à sympy. Mais une dernière chose importante à propos de l'import des fonctions de sympy.
Certains noms de fonctions de sympy, sin ou cos par exemple, sont identiques à ceux utilisés dans scipy. Pire, c'est vrai aussi pour des constantes comme pi ou e. Et il se trouve que les fonctions et variables sympy et scipy n'ont rien à voir ! Le pi de sympy est un symbole, le pi de scipy étiquette une valeur numérique de type float. Aussi, il faudra faire très attention chaque fois que vous utiliserez sympy avec d'autres modules comme scipy ou numpy, en utilisant correctement les namespaces, par exemple np.pi et sp.pi pour distinguer le pi de numpy et sympy respectivement.
Définir des variables et faire de l'arithmétique élémentaire
Définir les symboles
Imaginons que je veuille faire du calcul symbolique avec des fonctions dépendant de deux variables x et y. Sympy attend que je lui indique que x et y représentent des nombres, et même que je lui précise de quel type de nombre il s'agit : entier, rationnel, réel ou imaginaire. Sympy sait traiter, si vous lui dite, les nombres paris et impairs, les nombres premiers et l'infini.
Pour déclarer le symbole x réel, j'utilise l'instruction (je ne répète pas le prompt de la console sur les exemples à partir de maintenant):
x = Symbol("x")
Je n'ai rien spécifié concernant le type de la variable, le type "réel" étant le type par défaut. Si je veux déclarer un nombre imaginaire z, je taperai:
z = Symbol("z",imaginary=True)
Pour déclarer mes deux symboles x et y réels, je peux utiliser l'intruction suivante:
x,y = symbols("x,y")
Une remarque : il s'agit ici de réels au sens mathématique du terme, pas de flottants !
Pour l'exemple, déclarons un entier n :
n = Symbol("n",integer=True)
Evaluer, simplifier ou développer une expression
Au cours d'une séance de calcul symbolique, il est fréquent d'avoir à évaluer numériquement une expression littérale. Par exemple, au cours d'un calcul, j'ai besoin de connaitre la valeur numerique de \( \exp(pi)\). J'utilise la méthode evalf en lui indiquant le nombre de chifres que je veux, ici 10 :
exp(pi).evalf(10)
et j'apprends que cette valeur est 23,14069263. Vous constatez que 10 n'est pas le nombre de décimales, mais le nombre total des chiffres, avant et après la virgule.
Je peux aussi vouloir simplifier une expression. Prenons un exemple très connu. Soit r2 l'expression suivante:
r2 = sin(x)**2 + cos(x)**2
puis je demande à sympy de la simplifier en tapant:
simplify(r2)
il me répond la réponse attendue : 1 ! Cette fonction ressemble beaucoup à la fonction simplify de Maple ... Elle est assez riche et je vous invite à consulter le tutorial correspondant sur le site de sympy.
Autre exemple de simplification, celui de l'expression \( 2 \sin(x) \cos(x)\), avec l'instruction :
simplify(2*sin(x)*cos(x))
la réponse est \( \sin(2x) \). Assez utile lorsqu'on oublie ses tables trigo...
Vous pouvez aussi vouloir développer une expression, par exemple \( (x + 1)^2\). Je fais cela avec l'instruction :
expand((x +1 )**2)
j'obtiens bien la réponse attendue...
Un peu d'arithmétique
Quelques petites fonctions pas compliquées et très utiles.
Pour factoriser un polynôme :
factor(x**3 - x**2 + x -1)
qui me donne \( (x - 1)(x^2 + 1)\)
Réduire une expression au même dénominateur :
together(1/(1+x) + 1/y)
qui me renvoit \( \dfrac{x + y + 1}{y(x + 1)}\)
ou bien virer ces abominables racines carrées du dénominateur :
radsimp(1/(1 + sqrt(3)))
qui me retourne \( \dfrac{-1 + \sqrt(3)}{2}\)
Tracer la courbe représentative d'une fonction
Tracer une courbe 1D
Sympy utilise la librairie matplotlib de Python pour tracer des courbes, en bénéficiant de toutes les caractéristiques de cette librairie. Tracer une courbe est donc assez trivial. Voyons par exemple le tracé d'un sinus avec les axes et intervalles par défaut ([-10,10]). J'utilise la commande plot() :
plot(sin(x))
et j'obtiens :
Un autre exemple, le sinus cardinal \( \dfrac{\sin x}{x} \) entre \( -5 \pi \) et \( 5 \pi \):
plot(sin(x)/x, (x, -5*pi,5*pi))
que me donne :
Pas mal, vu l'économie de moyens mis en oeuvre...
Tracer une courbe 2D
Même principe, en utilisant la fonction plot3d(). Par exemple, tracer la surface \(x^2 + y^2 \) pour les intervalles (-5,5) pour x et (-5,5) pour y.
Il faut commencer par importer la fonction plot3d():
from sympy.plotting import plot3d
puis l'appeler avec les paramètres convenables:
plot3d(x**2+y**2, (x, -5,5), (y,-5,5))
que me donne :
Pour les différentes fioritures possibles, voir la doc de matplotlib !
Une curiosité : les zéros de la fonction zeta de Riemann
Il existe une librairie dans sympy qui s'appelle mpmath() et qui permet de tracer des fonctions complexes (au sens mathématique). Elle est un peu particulière, par exemple le nombre i n'existe pas et l'on ne peut pas employer le i de sympy. Il faut donc le déclarer comme ça :
I = complex(0,1)
Elle permet des choses géniales, par exemple tracer sur le plan complexe les racines nièmes de l'unité, par exemple ici la racine septième, avec l'instruction :
mpmath.cplot(lambda z:z**7-1,[-2,2],[-2,2], points =100000)
La fonction lambda me permet, en Python, de définir une fonction en ligne comme si j'écrivais def z : return z**7 - 1 . Attention , un peu long à afficher...
J'obtiens :
qui mérite aussi quelques commentaires:
- la gamme des couleurs dépend de la valeur de l'argument du nombre complexe. S'il s'agit d'un réel, la couleur est rouge.
- la transparence dépend du module du nombre complexe, plus c'est opaque et plus on tend vers 0 (noir pour 0)
Vous trouverez un peu plus d'éléments sur le site sympy déjà cité.
Enfin, pour le plaisir, le module mpmath contient une implémentation de la fonction zeta de Riemann (voir wikipedia pour ceux qui se demanderaient ce que c'est), qui permet de tracer les zéros de cette fonction, comme sur wiki ! Je commence par importer la fonction zeta():
from mpmath import zeta
puis je calcule et trace les zeros, attention le calcul est assez long :
mpmath.cplot(zeta,[-50,50],[-50,50], points = 100000)
et j'obtiens :
Je vous laisse consulter l'article de wikipedia pour l'interprétation de cette figure et les multiples secrets de la fonction zeta.
Sympy et les équations
Trouver les racines d'une fonction
C'est à dire résoudre l'équation f(x) = 0. Sympy fait ça très simplement. Je cherche par exemple les racines de \( f(x) = x^2 + x+ 1 \) :
roots(x**2 + x + 1)
sympy m'indique qu'il y a deux racines, l'une égale à \( -\dfrac{1}{2} - i \dfrac{ \sqrt 3}{2} \) et l'autre \( -\dfrac{1}{2} + i \dfrac{\sqrt 3}{2} \), ce qui est bien le résultat attendu.
Résoudre une équation
La fonction solve(), de sympy permet de résoudre les équations définies à l'aide de la fonction Eq(). Dans l'exemple ci-dessous, mon équation en x est \( x^2 = 2\). Je précise avec Eq() la valeur des deux termes de l'équation et j'indique à solve() la variable par rapport à laquelle je veux résoudre l'équation. Ici, ce n'est pas utile car il y a une seule variable, mais autant prendre les bonnes habitudes.
solve(Eq(x**2,2),x)
sympy retourne une liste Python qui contient les solutions, ici sans surprise \( -\sqrt 2\) et \( \sqrt 2\).
Résoudre un système d'équations
La résolution d'un système d'équations linéaires ou non est un jeu d'enfants! Il faut d'abord définir les équations. Dans mon exemple, il s'agit très classiquement de trouver x et y satisfaisant les conditions de somme et de différence :
equation1 = Eq(x + y, 50)
equation2 = Eq(x - y, 20)
puis j'utilise la fonction solve() pour chercher les solutions. Attention, les équations doivent être passées sous forme de liste Python à solve() :
solve([equation1, equation2])
le résulat est retourné sous forme de liste, ici x=35 et y=15 comme attendu.
Dériver et intégrer des fonctions
Chercher les limites d'une fonction
Trouver la limite d'une fonction avec sympy est encore un jeu d'enfant, avec la fonction limit(). Par exemple, cherchons la limite du sinus cardinal quand x tend vers 0. J'utilise l'instruction :
limit(sin(x)/x,x, 0)
et sans surprise, je trouve 1 ! Il suffit d'indiquer la fonction, la variable et la valeur vers laquelle on veut tendre...
On peut même cacluler la limite à droite et à gauche d'une fonction, par exemple pour 1/x, sa limite à droite:
limit(1/x,x,0,dir="+")
qui est égale à l'infini, et sa limite à gauche :
limit(1/x,x,0,dir="-")
égale à moins l'infini. Je suis sur que plein d'élèves de Première et de Terminale rêvent d'un tel outil :-)
Dériver une fonction
Fonction à une variable
Pour calculer une dérivée première, on utilise une fonction que nous avons déjà rencontré dans scipy, la fonction diff(). Par exemple, pour calculer la dérivée de sin(x) :
diff(sin(x),x)
qui nous retourne bien cos(x). Il est n'est pas indispensable ici de spécifier que l'on dérive selon la variable x, mais il est préfréable de la faire pour respecter l'ordre des argumesnt, ce qui sera utile ci-dessous. Un exemple un peu plus compliqué:
diff(sin(x)*exp(x)/x)
qui donne \( \dfrac{\exp x \sin x}{x} + \exp x \cos x - \dfrac{\exp x \sin x}{x^2} \). Si vous voulez vérifier...
Pour calculer une dérivée nième, rien de plus simple : il suffit d'indiquer l'ordre de dérivation (1 par défaut) à la suite de la variable :
diff(sin(x),x,2)
qui nous retourne -sin(x) comme attendu. On peut aussi utiliser l'écriture suivante, pour une dérivée seconde :
diff(sin(x),x,x)
Dérivée partielle
Soit la fonction \( f(x,y) = x^2 + y^2\), que je décris avec l'instruction suivante (lors de l'ouverture de session, le symbole f de fonction a été déclaré):
f = x**2 + y**2
Je calcule la dérivée partielle \( \dfrac{\partial f(x,y)}{x}\) avec :
diff(f,x)
qui me retourne 2x. Le principe est le même pour dériver par rapport à y ou pour obtenir une dérivée partielle nième.
On peut s'amuser à calculer la différentielle totale de f :
diff(f,x) + diff(f,y)
qui nous retourne \( 2x + 2y\).
Intégrer une fonction
Les intégrales simples
Le calcul d'intégrale indéfinie est aussi simple que le calcul de la dérivée d'une fonction. Soit par exemple la fonction \( g(x) = sin(x)*cos(x)\). Le symbole de la fonction g a été déclaré automatiquement par init_session(). Je définie la fonction g(x) et je calcule son intégrale indéfinie en utilisant l'instruction integrate() :
g = sin(x)*cos(x)
integrate(g,x)
sympy me retourne le résultat : \( \dfrac{1}{2} \sin^2 x \), que vous pourrez facilement vérifier manuellement en remarquant la nature de \( \cos x dx\)...
Il est tout aussi facile de calculer une intégrale entre deux bornes, par celle de sin(x) entre 0 et \( \dfrac{\pi}{2}\) :
g = sin(x)
integrate(g, (x,0,pi/2))
Notez la manière de définir les bornes !
Vous trouverez plus de détails sur l'intégration dans sympy sur le tutorial de sympy, en particulier sur les intégrales exponentielles, log et trigo, réelles ou complexes.
Les intégrales multiples
Reprenons notre fonction \( f(x,y) = x^2 + y^2\), pour calculer l'intégrale double selon x et y (calculer la surface d'un disque donc...). Je définie la fonction et calcule d'abord son intégrale indéfinie :
f = x**2 + y**2
integrate(f,x,y)
ce qui me donne \( \dfrac{x^3 y }{3} + \dfrac{x y^3 }{3}\)
Calculons miantenant cette intégrale pour x et y variant entre 0 et 1 :
integrate(f, (x,0,1),(y,0,1))
La méthode de déclaration des bornes du domaine d'intégration est la même que pour les intégrales simples.
Manipuler des matrices
Construire une matrice
On utilise l'instruction Matrix (c'est original) pour créer une matrice de dimension donnée et pour la remplir :
m1 = Matrix([[1,2,3],[4,5,6]])
J'ai construit ici la matrice m1, très simplement.
Opérations sur les matrices
Créons deux matrices carrées m1 et m2, afin de pouvoir faire quelques opérations. J'indique les intsrructions mais pas le résultat: vous le verrez bien et vous pourrez le vérifier manuellement, ce qui sera un excellent exercice !
m1 = Matrix([[0,1],[1,0]])
m2 = Matrix([[1,0],[0,1]])
une somme :
m1+m2
un produit par un scalaire
3*m1
un produit de matrices :
m1*m2
la transposée d'une matrice :
m1.transpose()
l'inversion d'une matrice
m1**-1
curieusement, il ne s'agit pas là d'une méthode, comme pour la transposée, mais d'un calcul... On peut d'ailleurs le généraliser à la puissance nième d'une matrice.
Mais il existe aussi une méthode, la méthode inv():
m1.inv()
Par défaut c'est la méthode du pivot de Gauss qui est utilisée, mais on peut utiliser la méthode LU avec :
m1.inv("LU")
le noyau d'une matrice
m1.nullspace()
Déterminant d'une matrice
Le calcul du déterminant d'une matrice :
m1.det()
Valeurs propres et vecteurs propres
Rappelons qu'un vecteur propre d'une matrice M est un vecteur v tel que \( M*v = \lambda*v\). Le scalaire \( \lambda \) est une valeur propre. Les valeurs propres d'une matrice sont les racines de son polynome caractéristique.
Considérons une matrice mp une matrice carrée 3x3:
mp = Matrix([[1,1,0],[0,1,0],[0,0,1]])
Son polynôme caractéristique peut être obtenu par la méthode charpoly() associée à la méthode as_expr() qui permet d'en donner l'expression symbolique :
mp.charpoly(x).as_expr()
Pour notre matrice, son polynôme caractéristique est \( x^3 - 3x^2 + 3x - 1\). On peut chercher ses racines :
roots(x**3 -3*x**2 + 3*x -1)
égales à 1 et 3.
Je vais maintenant calculer les valeurs propres de mp par la méthode eigenvals() de sympy avec :
mp.eigenvals()
qui me retourne 1 et 3. Nous avons bien vérifié, qu'au moins sur cet exemple, les valeurs propres d'une matrice sont les racines de son polynôme caractéristiques. Les vecteurs propres associé à ces deux valeurs propres sont obtenus par la méthode :
mp.eigenvects()
A titre d'exercice, vous pourrez vérifier avec sympy que les deux vecteurs propres [1,0,0] et [0,0,1], des vecteurs colonnes attention, satisfont bien à la définition. Il faut d'abord construire un vecteur colonne v avec les valeurs retournées par eigenvects() :
v = Matrix([1,0,0])
puis calculer le produit mp*v
mp*v
et v par la première valeur propre, ici 1
1*v
et comparer....
Résoudre un système linéaire donné sous forme matricielle
Imaginons un système linéaire de la forme A*X = b, avec A la matrice des coefficients et b le vecteur colonne des paramètres. Nous savons que X = Inverse(A)*B. Si nous reprenons notre système linéaire de tout à l'heure, nous avons les deux matrices A et B suivantes :
A = Matrix([[1,1],[1,-1]])
B = Matrix([50,20])
calculons le produit :
A.inv()*B
Nous obtenons bien le vecteur colonne X = [35,15]. Il est également possible d'utiliser l'instruction linsolve() comme ceci :
linsolve(Matrix(([1,1,50],[1,-1,20])),(x,y))
qui donne bien sur le même résultat.
Retour sur les opérateurs différentiels
Dans une autre page de TangenteX, nous abordons le calcul numérique de différentes opérateurs différentiels comme le gradient, la divergence ou le rotationnel. Sympy nous permet d'envisager leur calcul littéral, ce qui peut être fort utile pour certains exercices. Voyons comment procéder
Avant de commencer à travailler sur les opérateurs différentiels, vous devez importer la librairie vectorielle avec l'instruction :
from sympy.vector import *
Gradient
Nous allons débuter ce calcul par la définition d'un système de coordonnées cartésiennes, que nous appellerons R :
R = CoordSysCartesian('R')
Dans notre référentiel R, les vecteurs de base sont désignés par R.i, R.j et R.k. Un vecteur V (1,2) s'écrira par exemple V = R.i + 2*R.j. Dans le cas d'un champ, les grandeurs x, y et z sont désignées par R.x, R.y et R.z. Par exemple, un champ scalaire E d'équation \( E(x,y) = \dfrac{1}{x^2 + y^2} \), qui vous rappelle sans doute quelque chose, s'écrira sous sympy:
E = 1/(R.x**2 + R.y**2)
Cherchons son gradient à l'aide de l'instruction :
R.delop.gradient(E).doit()
Dans cette instruction, la méthode delop donne accès aux coordonnées du vecteur gradient et la méthode doit() évalue le résulat. Nous obtenons l'expression du gradient : \( \vec{grad} \: E = -\dfrac{2x}{(x^2 + y^2)^2} \vec{e}_i - \dfrac{2y}{(x^2 + y^2)^2} \vec{e}_j\).
La méthode est la même pour un champ vectoriel, mais dans sa définition, il faut prendre garde à bien indiquer les vecteurs de base, par exemple:
F = (R.x**2/2)*R.i - R.x*R.y*R.j
Divergence
Calculons la divergence du champ vectoriel \( F \) défini ci-dessus, la méthode dot() calculant la divergence :
R.delop.dot(F).doit()
Le résultat est celui attendu : la divergence de F est nulle. Si maintenant, nous calculons de la même façon la divergence du champ:
F = (R.x**2/2)*R.i + R.x*R.y*R.j
nous obtenons une divergence égale à \( div(\vec{F}) = 2x \), résultat également attendu !
Rotationnel
Soit le champ \( \vec{E} = -\sin (y) \vec{e}_i + \cos (x) \vec{e}_j \), dont le rotationnel est non nul. Définissons le champ sous sympy:
E = -sin(R.y)*R.i + cos(R.x)*R.j
puis calculons son rotationnel (curl en anglais)
curl(E,R)
Nous obtenons \( \vec{rot} \vec{E} = (\sin x + \cos y) \vec{e}_k \), un vecteur perpendiculaire au plan Oxy : normal le rotationnel est un produit vectoriel !
Contact : PhysiqueX ou 