Résolution de systèmes linéaires
| graphiboc |
Posté le 30 Sep 2009 à 13:32
|
|
 Messages : 262 GCPoints : 40352 |
Salut à tous, tout d'abord, je voudrais justifier la présence de mon post dans cette section : ce que je vais présenter ici ne sera pas une librairie proprement dite, mais un code qui pourrait par la suite en devenir une. L'objectif : Le problème ici s'intéresse à la résolution de système linéaire (ie du type 3x+5y-7z = 13) et plus particulièrement aux systèmes linéaires de grande taille. L'objectif c'est d'obtenir un code très performant sur ce type de système, en choisissant la méthode la plus appropriée pour la résolution. Le topic va donc consister tout d'abord à coder chaque méthode de résolutions, à en expliquer vaguement le principe (je ne ferai pas un cours de maths) et à enfin les rassembler en un seul bloc qui choisira automatiquement la méthode la plus appropriée. On s'interressera en particulier aux matrices creuses (qui ont beaucoup de 0) L'utilité : Sans doute que pour la plupart vous ne voyez pas d'interet dans la résolution de grands systèmes linéaires. Mais si vous vous intéressez un jour à des simulations physique non simplistes, vous en aurez peut etre besoin. En fait de tels systèmes peuvent permettre d'approcher la solution d'une équation différentielle assez complexe via une méthode de discrétisation que l'on appelle méthode des différences finies. En gros, vous découpez votre espace en un maillage de points, vous numérotez chaque point et vous cherchez à avoir une valeur approchée de votre solution en chacun de ces points, à partir d'une condition initiale qui est la connaissance de la valeur que prend la solution sur un contour de points. Plus il y a de points, plus c'est précis, mais plus le système est gros (ex : un maillage 10x10 fourni deja une centaine d'équations) Les méthodes : Elles seront regroupées en deux grandes catégories : les méthodes exactes et les méthodes itératives. Les méthodes exactes permettent d'obtenir la solution exacte du système. En revanche se sont les plus couteuses en terme de temps d'exécution. Les méthodes itératives permettent elles de converger vers la solution avec un écart maitrisé (donc d'obtenir une valeur approchée). De manière générale elle sont simples à effectuer, donc peu couteuses. Mais leur vitesse de convergence peut être lente. Ce qu'il faut bien comprendre c'est que ces méthodes ne peuvent pas forcément s'appliquer à n'importe quel système, d'où l'utilité d'analyser le système avant résolution, mais aussi que la méthode à utiliser pour un résultat optimal dépend directement du système. Une dernière remarque : a priori je vais coder tout cela sous Dark Basic (peut etre ensuite matlab ou C++), et donc même les méthodes dites exactes pourront donner des valeurs approchées voire légérement erronées. En revanche, le même algorithme utilisé dans un logiciel de calcul formel (ex : Mapple) donnera la solution exacte. La transposition d'un langage à l'autre ne devrait pas être trop difficile. Les méthodes exactes : -Cramer : la plus lourde, jamais utilisée, complexité en factorielle n (n est le nombre d'équation du système). -Pivot de Gauss : assez simple, mais complexité en n^3 (Je m'adresses ici en particulier à Darktib : tu découvriras cette méthode normalement cette année, et celle de Cramer l'année prochaine. Cramer est pratique pour les systèmes de petite taille, et les formules sont mathématiquement plus simple, mais difficiles à calculer...) -Factorisation LU : celle-ci est assez particulière. Son but, c'est de résoudre le même système, avec juste les seconds membres qui changent. L'interet c'est par exemple d'avoir un ensemble d'équations dont le second membre dépend du temps, on résout une première fis le système avec une complexité en n^3 puis on peut faire varier indéfiniment le second membre et résoudre le système en complexité n^2, ce qui est loin d'être négligeable. -Amélioration : les deux dernières méthodes peuvent être améliorées si la matrice du système est symétrique, symétrique définie positive, ou à bande. (on peut passer en complexité linéaire pour une matrice à bande, dans certains cas) Les méthodes approchées : -Jacobi -Gauss-Seidel Sur lesquels je ne m'étend pas pour le moment Matrice de système ? Peut être que certainsne voient pas ce que je veux dire par là. Je m'explique : Supposons que l'on dispose du système constitué des 3 équation suivantes : x+2y+5z = 3 y+3z = 1 x+z = 2 les seconds membres sont pour chaque équation, dans l'ordre : 3, 1, 2 les inconnues sont x, y, z la matrice du système est : 1 2 5 0 1 3 1 0 1 (dsl je n'ai pas pu expliquer plus pour des histoires de présentation) Actuellement je n'ai encore rien programmé. Le code arrivera donc au fur et a mesure. Malheureusement je ne pourrais pas expliquer vraiment le fonctionnement de ces méthodes, ca ne servirait qu'à écrire des formules inintéressantes. Mais les codes seront suffisamment simples pour que vous puissiez les reprendre dans le langage de votre choix. Les deux premières méthode (Gauss et Lu car Cramer ne présente aucun interet) seront bientot disponibles.
L'Homme est un créateur...
|
|
| Darktib |
Posté le 30 Sep 2009 à 19:05
|
|
 Messages : 4017 GCPoints : 347288 |
Ca m'a l'air très interessant et utile (ca peut aussi servir pour autre chose que des simulations physiques). Pivot de Gauss : vu il n'y a meme pas une semaine  et pour l'instant ca ne me sert toujours à rien: les systèmes que je résous sont simplissimes...D'ailleurs ca me fait penser qu'il faudra que tu ajoute une méthode pour des systèmes très simples, genre: L1 : 2x - y + 15/8 = 0 L2 : 3x + y = pi/2 (au pif), il suffit de faire L1 + L2 pour trouver x et d'utiliser une des deux lignes pour y. Donc ici, pas besoin d'algos complexes, et ca doit etre beaucoup plus rapide... Autre question : en quel langage ? C, C++, dbp, OCaml ? |
|
| graphiboc |
Posté le 01 Oct 2009 à 15:39
|
|
|
Messages : 262 GCPoints : 40352 |
Nous on l'avait fait un peu plus tard dans l'année je pense. Ca te servira pour déterminer facilement l'inverse d'une matrice, ou pour trouver des équation d'espaces vectoriels, etc. En fait l'algorithme du pivot de gauss consiste justement en ce genre de simplifications. Pour le langage j'ai adopté le dbp pr le moment. Lorsque tout sera fini, je ferai surement une transcription en C++. Mais ce sont des algorithmes très mathématiques, donc facilement portables d'un langages à l'autre. Cas particulier : matrice triangulaire supérieur Voici le premier code de résolution, qui traite d'un cas très particulier : le cas d'une matrice de système triangulaire supérieure. Typiquement une telle matrice correspond à un système de ce type : 3x + 2y + 8z = 4 5y + z = 1 z = 7 la matrice du système est : 3 2 8 0 5 1 0 0 1 je pense qu'ainsi vous voyez mieux ce que triangulaire supérieure signifie. Ce genre de système se résoud en cascade : on part de z = 7, puis de 5y + z =1 ce qui donne y puis la première équation donne x. On "remonte" le système. Ce cas particulier est très utile car les algorithmes de gauss et de LU visent à se ramener à ce genre de système (en fait c'est exactement ce qu'ils font). Voilà le code : Code :
global n : n=3 `taille du système
DIM A(n,n) as float `matrice du système
DIM X(n) as float `solution
DIM B(n) as float `second membre
`cas d'une matrice trigonale supérieure
function SolveUp()
local S as float
X(n) = B(n) / A(n,n)
for k = n-1 to 1 step -1
S = B(k)
for j = k+1 to n
S = S - A(k,j)*X(j)
next j
X(k) = S / A(k,k)
next i
endfunction
Voilà pour le moment. A noter que le code va peu à peu se remplir des différentes fonctions, je ne posterai donc que la nouvelle fonction à chaque fois, mais il s'agit d'un tout.
Dernière édition le 01 Oct 2009 à 15:54
L'Homme est un créateur...
|
|
| graphiboc |
Posté le 01 Oct 2009 à 15:47
|
|
|
Messages : 262 GCPoints : 40352 |
Cas particulier n°2 : matrice triangulaire inférieure De la même manière, on peut avoir une matrice de système triangulaire inférieure. par exemple le système suivant : x = 7 x + 2y = 2 5x + 3y + 2z = 3 donne pour matrice de système : 1 0 0 1 2 0 5 3 2 cette fois, on "descend" le système pour le résoudre. On peut donc écrire de même : Code :
`cas d'une matrice trigonale inférieure
function SolveDown()
local S as float
X(1) = B(1) / A(1,1)
for k = 2 to n
S = B(k)
for j = 1 to k-1
S = S - A(k,j)*X(j)
next j
X(k) = S / A(k,k)
next i
endfunction
Ce cas ci et le précédent seront utilisés dans la factorisation LU. Remarque : On peut voir que la notion de système à matrice triangulaire dépend du choxi de l'ordre des inconnues et des équations. Par exemple, dans le système précédent, en posant pour inconnues x'=y, y'=z et z'=x on obtient, en prenant l'ordre x',y',z', la matrice de système : 0 0 1 2 0 1 3 2 5 matrice qui n'est pas triangulaire. PS: tu as une idée d'autre application?
Dernière édition le 01 Oct 2009 à 15:54
L'Homme est un créateur...
|
|
| Darktib |
Posté le 03 Oct 2009 à 17:50
|
|
|
Messages : 4017 GCPoints : 347288 |
Je ne vois pas vraiment ce que tu veux dire avec 'application' - surtout que le terme à un sens mathématique (mais ca tu le sait^^). | |
| Mod |
Posté le 04 Oct 2009 à 13:28
|
|
 Messages : 4954 GCPoints : 2100823 |
Tout ça fait partiellement partie de mes connaissances (du moins en théorie, les cours de maths valent ce qu'ils valent), mais pour ma part je n'ai absolument aucune idée de ce à quoi peut servir cela en pratique. Je serais d'ailleurs bien curieux d'en voir une démonstration en pratique (même du point de vue d'un moteur physique, j'ai du mal à cerner où exploiter cela). |
|
| Darktib |
Posté le 04 Oct 2009 à 14:40
|
|
|
Messages : 4017 GCPoints : 347288 |
Raycasting : tu possède l'équation de plusieurs forme, et l'équation d'une droite. Tu veux déterminer le point d'intersection du rayon et de la forme. Dans ce cas, il faut résoudre un système linéaire. | |
| graphiboc |
Posté le 04 Oct 2009 à 15:24
|
|
|
Messages : 262 GCPoints : 40352 |
par "application" je voulais juste signifier "utilisation". Il s'agit ici de résoudre des systèmes à plusieurs centaines d'inconnus, donc le raycasting, oui et non. En fait, souvent en physique on dispose d'une équation différentielle, mais dont on ne connait pas la solution mathématique générale (équation de poisson par exemple). L'intérêt ici c'est d'approcher les valeurs de la solution sans résoudre l'équation différentielle. En gros, c'est une résolution numérique et non analytique. Le principe est de discrétiser les équations différentielles en utilisant le principe de base suivant : f'(x) = limite lorsque h tend vers 0 de [f(x+h) - f(x)]/h, où h est le pas du système. On discrétise en faisant sauter la limite et en posant directement f'(x) = [f(x+h)-f(x)]/h Plus h est petit, plus le résultat est précis. Pour calculer par exemple le champ magnétique dans un espace 2-dimensionnel, on découpe cette espace en maillage, en introduisant un certain nombre de noeuds (une espèce de quadrille carré de coté égal à h). En appliquant la formule de f' en chacun des points (à noter que l'on a cette fois-ci des dérivées partielles en x et en y), on obtient un système linéaire à résoudre. En numérotant correctement les noeuds et les équations, on peut obtenir une matrice tridiagonale par bloc, assez facile à résoudre. Dites moi si je n'ai pas été assez clair. Bon passons à la suite : Méthode de gauss : Je ne vais pas m'étendre sur la théorie qu'il y a derrière cette méthode (théorie assez simple cependant). Le but c'est de se ramener à un système triangulaire, que l'on pourra résoudre en cascade. Le principe est en fait assez simple : c'est exactement le même que lorsque vous avez une équation toute bête à deux inconnues x et y, et que vous exprimer x en fonction de y pour pouvoir injecter le résultat dans la 2e équation et obtenir ainsi une équation simple à une seule inconnue. En généralisant à n inconnues, on obtient un certain nombre de formules et c'est exactement ce que fait l'algorithme (à noter que cette méthode différe légérement de la méthode du pivot généralement enseignée). voici le code : Code :
`méthode du pivot de Gauss, sans permutation
function SolveGauss()
local r as float
for k=1 to n-1
for i=k+1 to n
r = A(i,k) / A(k,k)
for j=k+1 to n
A(i,j) = A(i,j) - r * A(k,j)
next j
B(i) = B(i) - r * B(k)
next i
next k
SolveUp()
endfunction
Cependant cette méthode n'est pas assez générale, d'où la précision "méthode sans permutation". En effet elle ne s'applique qu'à certaines matrice pour lesquelles on ne tombera jamais sur un 0 dans la diagonale au cours de l'algorithme. Pour remédier à cela, lorsque l'on tombe sur un zéro, il faut permuter avec une ligne inférieure de la matrice et continuer le processus (ceci est toujours possible si le système admet une solution). Ceci engendre pas mal de modification dans le code. Tout d'abord l'ajout d'un tableau de permutations. Le début du code devient : Code : global n : n=4 `taille du système global DIM A(n,n) as float `matrice du système global DIM X(n) as float `solution global DIM B(n) as float `second membre global DIM P(n) as integer `permutations Ce qui induit l'obligation de réinitialiser ce tableau avant de résoudre chaque système : Code : `initialiser les données nécessaires à la résolution (à utiliser si les équations changent) function Init() `réinitialiser le tableau de permutation for i=1 to n P(i)=i next i endfunction Ne pas oublier une fonction qui permet de permuter 2 lignes : Code : `permute les lignes i et j function Permuter(i,j) tmp=P(i) P(i)=P(j) P(j)=tmp endfunction et enfin il faut modifier la fonction en prenant en compte les permutations : Code :
`méthode du pivot de Gauss, avec permutation
function SolveGaussPerm()
local r as float
PrintMatrice()
for k2=1 to n-1
k=P(k2) `prise en compte de la permutation
`Vérification de permutation
if A(k,k2)=0
`on cherche une ligne afin d'avoir ce coefficient non nul
for l=k+1 to n
if A(l,k2)<>0
Permuter(k,l) `on permute
exit `on passe à la suite
endif
next l
endif
for i2=k2+1 to n
i=P(i2) `prise en compte de la permutation
r = A(i,k2) / A(k,k2) `pas de permutation sur les colonnes, d'où le k2!
for j=k2+1 to n
A(i,j) = A(i,j) - r * A(k,j)
next j
B(i) = B(i) - r * B(k)
next i
PrintMatrice()
next k
SolveUp()
endfunction
Désolé de ne pouvoir m'étaler plus longtemps sur la théorie, essayez de comprendre le code en admettant la partie mathématique. Un petit exemple : On veut résoudre le système suivant x+z+t=8 x+2y+3z+4t=30 x+3y+6z+10t=65 x+4y+10z+20t=119 La résolution mathématique donne x=1,y=2,z=3,t=4 le code pour résoudre ce système est le suivant : Code : global n : n=4 `taille du système global DIM A(n,n) as float `matrice du système global DIM X(n) as float `solution global DIM B(n) as float `second membre global DIM P(n) as integer `permutations A(1,1)=1 : A(1,2)=0 : A(1,3)=1 : A(1,4)=1 A(2,1)=1 : A(2,2)=2 : A(2,3)=3 : A(2,4)=4 A(3,1)=1 : A(3,2)=3 : A(3,3)=6 : A(3,4)=10 A(4,1)=1 : A(4,2)=4 : A(4,3)=10 : A(4,4)=20 B(1)=8 : B(2)=30 : B(3)=65 : B(4)=119 Init() SolveGaussPerm() print X(1) print X(2) print X(3) print X(4) sync wait key end on obtient bien x=1,y=2,z=3,t=4
Dernière édition le 04 Oct 2009 à 15:35
L'Homme est un créateur...
|
|
Voir le profil
Envoyer un message privé