Découvrir Scilab/Résolution d'équations
Résolution d'équations
Scilab dispose de plusieurs primitives permettant la résolution d'équations. Nous présentons ci-après quatre fonctions internes, mais il en existe d'autres pour des systèmes spécifiques.
Équation polynomiale.
Rappelons que la commande codice_1 détermine les racines du polynôme "p", c'est-à-dire que cette fonction résout l'équation
L'équation
se résout donc avec
solution = roots(p - y)
Équation linéaire matricielle.
La commande codice_2 permet de déterminer le noyau de la matrice, c'est-à-dire, pour une matrice M, de résoudre l'équation matricielle
la syntaxe est simplement :
X = kernel(M)
La division de matrice consiste à résoudre une équation linéaire :
Notons que l'on peut représenter un système d'équations linéaires par une matrice ; la division matricielle est donc un moyen de résoudre un tel système.
Si la matrice A est creuse, il est recommandé d'utiliser codice_6 ; pour résoudre, codice_4 :
X = umfpack(A, "\" , B)
Par ailleurs, les résolutions numériques d'équations matricielles font souvent intervenir des décompositions d'une matrice :
Système d'équations linéaires.
Scilab permet la résolution numérique d'un système d'équations linéaires, avec la fonction codice_11
[x, kerA] = linsolve(A, b)
où A est la matrice réelle des coefficients du système d'équations, et "b" un vecteur de constantes. Les solutions de
sont de la forme
"w" étant un réel arbitraire.
Exemple
-->A = [3 1 ; 4 -1];
-->b = [-5 ; -9];
-->linsolve(A, b)
 ans =
 2.
 - 1.
-->[x0, kerA] = linsolve(A,b)
 kerA = 
 x0 =
 2.
 - 1.
Exemple
-->A = [3 1] ; b = -5;
-->linsolve(A, b)
 ans =
 1.5
 0.5
-->[x0, kerA] = linsolve(A,b)
 kerA =
 - 0.3162278
 0.9486833
 x0 =
 1.5
 0.5
Si le système n'admet pas de solution, Scilab affiche le message d'erreur : codice_12.
Scilab permet également la résolution symbolique d'un système linéaire avec la fonction codice_13
w = solve(A, c)
où A est une matrice triangulaire supérieure de réels (les coefficients du système d'équation), "c" est un vecteur de chaînes de caractères (la partie à droite des signes égal), et "w" est la matrice résultat de
Exemple
-->A = ["1", "1" ; "0", "1"] ; c = ["a1" ; "a2"];
-->solve(A,c)
 ans =
!a1-a2 !
!a2 !
Scilab propose aussi la méthode du gradient conjugué avec la fonction codice_14. La syntaxe élémentaire est :
X = conjgrad(A, b)
Système d'équations quelconque.
Scilab permet la résolution numérique d'un système d'équations quelconques, avec la fonction codice_15.
Considérons un système de "n" équations à "n" inconnues ("x"1, …, "xn"), de la forme
On définit, en tant que fonction externe (voir "Programmation > Fonctions") la fonction ƒ qui associe à un vecteur x("x"1, …, "xn") un vecteur v("vx"1, …, "vn") :
Le système d'équations consiste donc à résoudre : ƒ(x) = 0.
La syntaxe utilisée est :
[x] = fsolve(x0, f)
où codice_16 est la fonction externe telle que définie ci-dessus et codice_17 une estimation de la solution. Dans le cas d'un "système" d'équations, l'entrée x est donc un vecteur de dimension "n" et la fonction codice_16 renvoie un vecteur de dimension "n".
Scilab utilise un algorithme du type méthode hybride de Powell. La fonction peut donner le vecteur v correspondant à la solution x estimée (v devant être proche de 0), avec la syntaxe :
[x, v] = fsolve(x0, f)
Il est recommandé de fournir la jacobienne de ƒ, sous la forme d'une fonction externe ƒJ :
où J est la matrice jacobienne :
[x] = fsolve(x0, f, fJ)
Exemple
function [y] = f(x)
 y = [3, 1 ; 4, -1]*x - [5 ; 9];
endfunction
x0 = [0 ; 0];
solution = fsolve(x0, f)
Si la fonction est paramétrée, alors
Exemple
function [Y] = f(X, A, b)
 Y = A*X + b;
endfunction
param1 = [1, 0 ; 0, 3];
param2 = [2 ; 4];
X0 = [0 ; 0]
solution = fsolve(X0, list(f, A = param1, b = param2));
// alternative : solution = fsolve(X0, list(f, A = param1, b = param2));
disp(solution);
Exemple complet.
Considérons un système mécanique simple : une poutre pivote autour d'un point fixe A("x"A, "y"A) ; un vérin fixe pousse sur cette poutre selon un contact ponctuel au point B("x"B, "y"B). La hauteur de la poutre vaut "e" = , le pivot est implanté à mi-hauteur. Le vérin est implanté à une distance horizontale L = . On a donc "x"B – "x"A = L, et l'altitude "y"B varie. Le sujet présente une symétrie plane, on se place donc dans le plan de symétrie (A, "x", "y"). Dans ce repère, on a donc :
On désire connaître l'inclinaison de la poutre en fonction du déplacement du vérin. Pour cela, on détermine l'équation de la droite génératrice de la poutre, située sur sa face inférieure. L'équation générale de cette droite (D) est :
et l'on veut donc déterminer les paramètres "a" et "b". Les conditions géométriques sont :
On doit donc résoudre un système de deux équations, qui n'est pas linéaire.
On désire développer du code que l'on pourra copier-coller. Donc, au niveau des notations :
Le système d'équations devient donc :
Or ce problème présente deux solutions, puisqu'il existe deux droites passant par B et tangentes au cercle de centre A et de rayon "r". Il faut donc :
Pour accélérer les calculs, lorsqu'une expression revient plusieurs fois, on la met dans une variable ; on rappelle que seul "y"B = A(4) est variable. On identifie donc :
soit
Donc,
On a donc deux vecteurs de valeurs "a" et "b", les valeurs "a"("i") et "b"("i") correspondant au déplacement du piston "y"B("i"), ce qui permet de déterminer l'inclinaison de la poutre pour chaque valeur "y"B("i"). La valeur "b"("i") ne nous sert à rien ("b" est juste un intermédiaire de calcul permettant de résoudre le système), mais nous la conservons tout de même car elle permet de tracer rapidement la droite solution si cela nous intéresse.
Une fois ceci fait, on désire déterminer les effort aux liaisons, c'est-à-dire la force FA que l'axe exerce sur la poutre en A, et la force FB que le vérin exerce sur la poutre en B. La force FA est totalement inconnue :
et on connaît la direction de la force FB : c'est le vecteur normal (–"a", 1). Le vecteur peut donc s'écrire sous la forme :
Le système a une "m", le centre des masses étant en G("x"G, "y"G) ; prenons arbitrairement "m" = , "x"G = et "y"G = . L poids P s'écrit p(0, –"mg") avec "g" = .
Les équations de la statique s'écrivent
soit
Si l'on définit le vecteur des inconnues
et les vecteurs de paramètre
les équations de la statique s'écrivent sous la forme d'une équation matricielle :
En résolvant ce système, on peut en déduire les intensités des forces :
Concrètement, le programme développé est :
// paramètres du problème géométrique
xA = 0;
yA = 0;
e = 50; // en mm
L = 1000; // en mm
xB = L;
course = 200; // course du vérin en mm
yBmin = -e/2;
yBmax = yBmin + course;
// paramètres du problème de statique
xG = 700; // en mm
yG = 200; // en mm
m = 100; // en kg
g = 9.81; // en m/s^2
n = 11; // nombre de points calculés
// paramètres dérivés
r = e/2;
yB = linspace(yBmin, yBmax, n); // vecteur de valeurs de y pour B
// fonction
function [Y] = contraintes_geo(X, A)
 // droite d'équation y + az + b = 0
 // a et b à déterminer
 // X(1) = a ; X(2) = b
 // A(1) = xA ; A(2) = yA ; A(3) = xB ; A(4) = yB ; A(5) = r ;
 // 1. B est sur la droite
 // yB + axB + b = 0
 Y(1) = A(4) + A(3)*X(1) + X(2);
 // 2. La distance de A à la droite vaut r
 // (yA + axA + b)^2/(1 + a^2) - r^2 = 0
 Y(2) = (A(2) + A(1)*X(1) + X(2))^2/(1 + X(1)^2) - A(5)^2;
endfunction
// * Programme principal *
// ********** Problème géométrique **********
// initialisation
param = [xA ; yA ; xB ; yBmin ; r]; // correspond à A dans la fonction
a = zeros(n);
b = a;
// simplification du calcul
facteur1 = 1/(param(3) - param(1));
facteur2 = (param(5) - param(2))*facteur1;
// recherche des valeurs
for i = 1:n
 param(4) = yB(i);
 facteur3 = param(4)*facteur1 + facteur2;
 X0 = facteur3*[-1 ; param(1)] + [0 ; param(2)]; // point de départ
 solution = fsolve(X0, list(contraintes_geo, param));
 pente = -solution(1);
 a(i) = -pente;
 b(i) = solution(2);
 // vérification que l'on a pris la bonne solution
 penteAB = (param(4) - param(2))*facteur1;
 if pente < penteAB then // si c'est la mauvaise solution
 a(i) = %nan;
 b(i) = %nan;
 end
end
// affichage de l'angle
angle = atan(-a); // en radians
disp("position (mm) ; angle (°)");
disp([yB', 0.1*round((1800/%pi)*angle)]); // angle en °, arrondi
// ********** Problème de statique **********
// initialisation
FA = zeros(n, 2); 
C = zeros(n, 1);
iFA = C; // intensités des forces
iFB = C;
AF = [1, 0, -a(1)
0 1 1
0 0 L];
bF = -m*g*[0 ; 1 ; xG];
// Résolution des équations de la statique
for i=1:n
 AF(1, 3) = a(i);
 XF = linsolve(AF, bF);
 FA(i, :) = XF(1:2);
 C(i) = XF(3);
 iFA(i) = norm(XF(1:2));
 iFB(i) = C(i)*norm([a(i) ; 1]);
end
// affichage des intensités
disp("position (mm) ; FA (N) ; FB(N)");
disp([yB', round(iFA), round(iFB)]); // angle en °, arrondi
On obtient donc :
position (mm), angle (°)
 -25. 0. 
 -5. 1.1
 15. 2.3
 35. 3.4 
 55. 4.6
 75. 5.7
 95. 6.9
 115. 8.
 135. 9.1
 155. 10.2
 175. 11.3
 position (mm) ; FA (N) ; FB(N)
 -25. 294. 687.
 -5. 295. 687.
 15. 296. 687.
 35. 297. 688.
 55. 299. 689.
 75. 302. 690.
 95. 306. 692.
 115. 310. 693.
 135. 314. 695.
 155. 319. 698.
 175. 325. 700.
Notes.
"Calcul différentiel et intégral" < ↑ > "Optimisation d'une fonction"