Résoudre un système linéaire

Les systèmes linéaires en physique

La modélisation de nombreux systèmes physiques donne lieu à l'apparition d'équations linéaires de plusieurs inconnues, formant un système linéaire. Les premiers systèmes linéaires que l'on voit apparaitre en Première ou Terminale viennent en électrocinétique ou en statique. Puis viennent les problèmes d'électronique (ah les matrices de transfert des quadripôles!), d'hydraulique et de calcul des contraintes. Bref, dans presque toutes les branches de la physique, nous sommes amené à manipuler un système d'équations linéaires.
Rappelons d'abord à quoi ressemble un système linéaire (j'abrège...). C'est un ensemble d'équations que l'on peut écrire sous la forme:

a₁₁x₁ +a₁₂x₂ + a₁₃x₃ + ... + a_1nx_n = b₁
a₂₁x₁ +a₂₂x₂ + a₂₃x₃ + ... + a_2nx_n = b₂
................
a_n1x₁ +a_n2x₂ + a_n3x₃ + ... + a_nnx_n = b_n

Pour un système linéaire, on distingue trois caractéristiques:

le nombre d'équations n,
le nombre d'inconnues p,
le rang du système r, qui est égal au nombre d'équations indépendantes dans le système.

Le cas le plus favorable est celui où n=p=r, c'est à dire qu'il y a autant d'équations que d'inconnues et que toutes les équations sont indépendantes. Dans ce cas, la solution du système est unique et calculable simplement grâce au calcul matriciel. C'est le cas du système que j'ai écris ci-dessus.
Pour résoudre ce genre de système, on construit la matrice A formée des coefficients a_ij du système, le vecteur x formé des inconnues x_n et du vecteur b formé des seconds membres b_n, ce qui nous permet d'écrire notre système sous la forme matricielle Ax = b.
Le calcul matriciel nous apprend que si le déterminant de la matrice A est non nul, elle admet une matrice inverse notée A^-1, et que la solution du système est donnée par x = A^-1b.
La résolution d'un système linéaire de ce type revient donc à calculer une matrice inverse et à faire un produit d'une matrice par un vecteur. Si le rang du système est faible (2 ou 3 au maximum..), on peut faire les calculs à la main. Mais en physique, on se retrouve vite avec des systèmes de rang 10 ou 20, voire plusieurs centaines. En économie, on arrive à plusieurs milliers! Dans ce cas, il ne reste plus que l'ordinateur!
Si le déterminant est nul, c'est que le rang de la matrice (du système..) est inférieur à n, autrement dit que toutes les équations ne sont pas indépendantes. Et il faut arranger son système pour que cela devienne le cas, par exemple en fixant les valeurs des n-r inconnues. C'est souvent le cas en physique. Rappelons quelques règles, que l'on nous enseignait dans le temps en Spéciales:

si n > p, il y a un problème de modélisation, trop d'équations ! Il faut revoir l'application des lois utilisées pour la modélisation
si n < p, il manque des équations ! Soit vous avez oublié une ou plusieurs équations, soit le problème n'est pas soluble sauf à introduire un ou plusieurs paramètres qui permettront d'établir les équations manquantes.
si r < n, revoir la dépendance des équations! Il y a des équations écrites plusieurs fois...

Ce qui fait que l'établissement du système linéaire est souvent plus difficile que sa résolution...

Nous allons examiner dans la suite un système linéaire simple (et idéal avec n=r=p) et les différentes méthodes pour le résoudre.

Pour ceux que le sujet intéresse, je leur recommande la lecture du chapitre 4 du "Les maths en physique" de JP Provost et G.Vallée (université de Nice-Sophia Antipolis), qui détaille beaucoup d'exemples d'usage des systèmes linéaires en physique classique et quantique.

Le système linéaire pour l'exemple

Pour illustrer les différents programmes que je vous proposerai dans la suite, il me fallait trouver un système d'équations linéaires. Pas trop compliqué, connu de tous et de préférence linéaire!

Mon choix s'est porté sur un exemple ultra-classique et sujet bien connu d'exercices et de problèmes: le réseau de mailles. Voici donc un réseau de résistances alimentées par deux générateurs. Les fem V1 et V2, ainsi que les résistances R1 à R5 sont données. Le but du jeu consiste à calculer les intensités i1, i2 et i3 des courants circulant dans les différentes mailles.

Un célébrissime théorème bien connu de tous (lequel?) nous permet d'écrire les équations suivantes:
R1i1 + R2(i1 - i2) - V1 = 0
R2(i2 - i1) + R3i2 + R4(i2 - i3) = 0
R4(i3 - i2) + R5i3 + V2 = 0

Les intensités i1, i2 et i3 sont nos inconnues. Je vais donc réarranger ces trois équations pour en tenir compte, et je vais en profiter pour créer le second membre de mon système, ce qui me donne le système suivant:
(R1 + R2)i1 - R2i2 + 0i3 = V1
-R2i1 + (R2 + R3 + R4)i2 - R4i3 = 0
0i1 - R4i2 + (R4 + R5)i3 = -V2

Voilà donc le système d'équations à résoudre. Vous constatez que ces équations sont linéaires. Il est idéal: 3 équations indépendantes et 3 inconnues!
Sous forme matricielle, notre système s'écrit très simplement:
      A*C = V
en définissant la matrice A comme:
      (R1 + R2) -R2                   0
A = -R2            (R2 + R3+ R4)   -R4
       0             -R4                    (R4 + R5)
le vecteur C comme:
       i1
C = i2
       i3
et le vecteur V comme:
       V1
V = 0
       -V2
Nous obtenons donc un produit d'une matrice carrée (3x3) par un vecteur (3), qui est égale à un vecteur (3). Le but du jeu est bien sur de calculer la valeur des éléments i1, i2 et i3.

Dans la suite nous allons programmer plusieurs algorithmes qui permettent de faire ce calcul. Ces programmes sont très fortement inspirés des routines disponibles dans l'excellent ouvrage déjà maintes fois cité "Numerical Recipes in C", notre petit livre rouge... Il existe aussi pour le Fortran, mais il est bleu!

La méthode du pivot de Gauss

C'est sans doute la méthode la plus simple pour résoudre un système d'équations linéaires ! Au lycée, pour résoudre un système de 2 équations à 2 inconnues, on exprime dans la première équation, l'inconnue x1 en fonction de x2. Puis dans la deuxième équation, on substitue l'expression trouvée à la valeur x1, puis on calcule x2. Et enfin, connaisant x2, on calcule x1. Rien de plus simple...
Et bien, la méthode du pivot de Gauss est une généralisation de cette méthode ! Plus précisement, elle consiste à transformer la matrice A du système en une matrice triangulaire supérieure, qui permet de calculer x_n puis x_n-1 jusqu'à x₁. Voyons l'algorithme de façon un peu plus détaillée:

On commence par exprimer une inconnue x_i en fonction de toutes les autres en divisant par son coefficient a_1i. En principe, on commence par x₁, dont le coefficient a₁₁ doit bien sur être non nul. On appelle ce coefficient le pivot. Il faut cependant se méfier ! Le coefficient a₁₁ peut être très petit et donc provoquer l'apparition de valeurs très grandes dans les autres coefficients et par là des instabilités numériques. Il faut donc ne pas se précipiter sur le premier coefficient et choisir le coefficient pivot a_pi le plus grand. Celui-ci trouvé, on construit donc une matrice dans laquelle la ligne de l'inconnue x_p n'apparait plus.
on itère l'algorithme ci-dessus jusqu'à ce que la matrice soit triangulaire supérieure (voir cours de calcul matriciel pour ceux qui ne savent pas ce qu'est une matrice triangulaire supérieure...)
lorsqu'on a obtenu la matrice triangulaire supérieure, on calcule x_n puis les autres x_i de proche en proche.

Voilà ! L'algorithme est simple, mais pas très performant ni très stable numériquement. Le problème tient dans le choix du pivot, surtout pour les grands systèmes. De plus, il est assez gourmand en calculs. On va donc chercher autre chose...

Le code du programme et du projet Dev-Cpp associé sont téléchargeables. Le code est décomposé en trois grandes parties:

les routines utilitaires d'allocation et de libération de la mémoire utilisée par les matrices et les vecteurs. Ces routines sont ultra-classiques, utilisées dans beaucoup de programmes et très fortement inspirées de celles disponibles dans les "Numerical Recipes in C", qu'on ne présente plus!
la routine de calcul du pivot de Gauss et de résolution du système, l'élément original du programme. Elle implément strictement l'algorithme décrit ci-dessus.
le main, dont le rôle principal est de saisir les données et d'afficher les résultats...

L'utilisation du programme est très simple: vous le lancez en mode commande sous Windows (ou sous Linux), vous introduisez les données comme indiqué et le résultat s'affiche quasi instantanément.

La décomposition LU

Cette méthode est basée sur un constat de calcul matriciel: il est plus simple d'inverser une matrice triangulaire, qu'elle soit supérieure ou inférieure. Le principe est donc de décomposer le système Ax = b en deux systèmes plus simples à manipuler Ly = b et Ux = y. L est une matrice triangulaire inférieure (Low..) et U une matrice triangulaire supérieure (Up).

La première opération, la plus délicate, consiste à calculer les matrices L et U (routine FactorisationLU dans le programme DecompositionLU). Il existe un algorithme adéquat pour faire ça, c'est l'algorithme de Crout. Il réagence les lignes de la matrice A pour obtenir deux matrices triangulaires L et U. Je ne rentre pas dans le détails de cet algorithme ici. Le "Numerical recipes in C" donne un très bonne explication de cet algorithme au §2.3 de l'édition de 1989.

La seconde opération est plus simple: résoudre les deux systèmes Ly = b et Ux = y. La routine SolutionLU du programme fait ça :

on traite d'abord Ly = b pour calculer le vecteur y (on appelle cette étape de calcul la "forward substitution" dans les Recipes)
puis on traite Ux = y, qui nous donne le vecteur x (c'est la "backward substitution")

Le code du programme et du projet Dev-Cpp associé sont téléchargeables. Vous retrouverez la même structure que dans le programme PivotGauss. Seule diffère la routine de calcul, qui assure la décomposition LU. La saisie des données et l'affichage des résultats sont identiques.

Utiliser les logiciels de calcul

Résoudre un système linéaire avec Scilab

Scilab, construit sur une algèbre des matrices, est un outil particulièrement approprié pour traiter les systèmes d'équations linéaires. Voici un petit programme, dont le code est téléchargeable ici, qui résoud notre système:

//**************************************************************************

// Résolution du système linéaire A*C = V

// Dominique Lefebvre Décembre 2012

// TangenteX.com

//**************************************************************************

// Saisie des données

R = input("Valeur des resistances R1 à R5 [r1 r2 r3 r4 r5]: ");

U = input("FEM V1 et V2 [v1 v2]: ");

// Définition des matrices A et V

A = [R(1)+R(2), -R(2), 0; -R(2), R(2)+R(3)+R(4), -R(4); 0, -R(4), R(4)+R(5)];

V = [U(1); 0; -U(2)];

// Calcul de C

C = A\V;

// affichage des résultats

disp("Valeurs des courants i1, i2 et i3: ");

disp(C);

Le principe du programme est simple: on défini les matrices de données A et V, puis on calcule C.

Vous constatez que l'essentiel du programme est consacré à la saisie des données, la déclaration des matrices A et V et l'affichage des résultats. En fait, le calcul de C tient en une seule ligne, 5 caractères pour être précis! L'opérateur Scilab \ procède à la division matricielle à gauche. Dans le cas d'une matrice carrée régulière, ce qui est notre cas, cet opérateur équivaut à écrire C = inv(A)*V avec inv désignant la matrice inverse de A. Pour la petite histoire, Scilab utilise la librairie Lapack pour les matrices réelles et complexes.

Attention lors de l'usage du programme à entrer correctement les données! Par exemple, les valeurs de résistance doivent être entrées ainsi : [10 100 10 20 50]

Résoudre un système d'équations linéaires avec Scilab est donc un jeu d'enfant, du moins tant que le système n'est pas trop "tordu", mais ces cas ne rentrent pas dans notre scope...

Résoudre un système linéaire avec Maple

Voyons maintenant comment utiliser Maple pour traiter ce petit problème. Lancez Maple, puis lorsque vous aurez le prompt de Maple, entrez les commandes suivantes:

>R1 := 1:

>R2 := 1:

>R3 := 1:

>R1 := 1:

>R4 := 1:

>R5 := 1:

>V1 := 5:

>V2 := 5:

>fsolve({(R1+R2)*i1-R2*i2 = V1, -R2*i1+(R2+R3+R4)*i2-R4*i3 = 0, -R4*i2+(R4+R5)*i3 = -V2}, {i1, i2, i3});

{i1 = 2.500000000, i2 = 0., i3 = -2.500000000}

Les premières commandes ont pour objet l'initialisation des variables. Bien sur, vous choisirez vos valeurs pour les résistances et les FEM. j'ai mis ce qui me passait par la tête...
Puis vient la commande fsolve qui permet de résoudre numériquement un système d'équations. Dans cette commande, vous devez définir vos équations (ici 3 équations) puis vos variables (ici i1,i2 et i3). Vous appuyez sur ENTER et le tour est joué: Maple vous affiche la solution! Comparez donc avec les programmes ci-dessus pour les mêmes données en entrée, pas de problème!
Bien sur, rien ne vous interdit d'écrire un petit programme Maple pour faire ça, mais est-ce vraiment utile?

Résoudre un système linéaire avec Python

Le module NumPy (et aussi Scipy) de Python dispose d'une fonction de résolution des systèmes linéaires très performante et aussi simple de mise en oeuvre que Maple ou Scilab. Pour l'exemple, voici un code qui résoud notre système linéaire:

import numpy as np

# Définition des termes du système

R1 = 1.0

R2 = 1.0

R3 = 1.0

R4 = 1.0

R5 = 1.0

V1 = 5.0

V2 = 5.0

# Définition des matrices

A = np.matrix([[R1+R2,-R2,0],[-R2,R2+R3+R4,-R4],[0,-R4,R4+R5]])

B = np.array([V1,0,-V2])

# Résolution du système linéaire et affichage

print np.linalg.solve(A,B)

La réponse est bien sur identique à celles obtenues avec Scilab et Maple.

[ 2.5 0. -2.5]

Les moyens mis en oeuvre sont d'ailleurs quasiment identiques, puisqu'en python aussi on utilise une librairie matricielle et une méthode de décomposition en LU.

Pour conclure

Il existe bien d'autres algorithmes de résolution d'un système linéaire, plus efficaces et plus élégants. Vous les trouverez dans vos manuels de cours d'analyse numérique. Du point de vue de la programmation, le "Numerical Recipes in C" contient toutes les routines C pour résoudre ces systèmes. Attention cependant: dans les "Recipes in C", les indices de matrice et de vecteur vont de 1 à n, comme en Pascal et en Fortran. Pour ma part, je respecte les habitudes des programmeurs C en faisant partir mes indices de 0 à n-1.
Les programmes décrits ci-dessus ont une vocation essentiellement pédagogique. Ils ne sont pas très efficaces, voire carrément insuffisants dans les cas de grands systèmes ou de systèmes un peu "tordus". En cas d'utilisation sérieuse, je vous conseille vivement d'utiliser les routines des librairies de calcul bien connues comme LaPack, Root, Slatec ou GSL.
Je n'ai pas abordé ici les problèmes de conditionnement des systèmes, technique qui vise à manipuler la matrice d'un système afin d'améliorer la stabilité numérique de la solution. Là aussi, vous trouverez des informations sur le sujet dans votre manuel d'AN.