|
<noinclude>{{NavTitre|book={{BASEPAGENAME}}|prev=Calculs élémentaires|next=MatricesRésolution creusesd'équations}}</noinclude>
{{Scilab}}
Les commandes <code>conv</code> et <code>convol</code> n'utilisent pas le même algorithme. La fonction <code>convol</code> utilise la transformée de Fourier rapide, et est mieux adaptée lorsqu'il y a beaucoup de valeurs.
== Résolution d'équations ==
Scilab dispose de plusieurs primitives permettant la <font id="solveur1">résolution d'équations</span>. 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 <code>roots(p)</code> détermine les racines du polynôme ''p'', c'est-à-dire que cette fonction résout l'équation
: ''p''(''x'') = 0.
L'équation
: ''p''(''x'') = ''y''
se résout donc avec
<source lang="scilab">
solution = roots(p - y)
</source>
=== Équation linéaire matricielle ===
La commande <code>kernel()</code> permet de déterminer le [[w:fr:noyau (algèbre)|noyau]] de la matrice, c'est-à-dire, pour une matrice M, de résoudre l'équation matricielle
: MX = 0
la syntaxe est simplement :
<source lang="scilab">
X = kernel(M)
</source>
La division de matrice consiste à résoudre une équation linéaire :
* à droite : <code>X = A/B</code> est la solution de XB = A,
* à gauche : <code>X = A\B</code> est la solution de AX = B,
* on a <code>B/A = (A' \ B')'</code> ;
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 <code>umfpack()</code> ; pour résoudre, <code>X = A\B</code> :
<source lang="scilab">
X = umfpack(A, "\" , B)
</source>
Par ailleurs, les résolutions numériques d'équations matricielles font souvent intervenir des décompositions d'une matrice :
* <code>lu(M)</code> : factorisation LU d'une matrice carrée, décomposition en un produit d'une matrice triangulaire inférieure L et d'une matrice triangulaire supérieure U,
* <code>qr(M)</code> : factorisation QR, décomposition en un produit d'une matrice orthogonale Q et une matrice triangulaire supérieure R,
* <code>svd(M)</code> : décomposition en valeurs singulières ''(singular value decomposition)''.
=== Système d'équations linéaires ===
Scilab permet la résolution numérique d'un système d'équations linéaires, avec la fonction <code id="linsolve1">linsolve</code>
<source lang="scilab">
[x, kerA] = linsolve(A, b)
</source>
où A est la matrice réelle des coefficients du système d'équations, et ''b'' un vecteur de constantes. Les solutions de
: A×''x'' + ''b'' = 0
sont de la forme
: ''x''<sub>0</sub> + ''w''×kerA,
''w'' étant un réel arbitraire.
'''Exemple'''
: Le système d'équation
:: <math>\begin{cases}
3x+y=5\\
4x-y=9
\end{cases}
</math>
: est décrit par
:: <math>\mathrm{A} = \begin{pmatrix}
3 & 1 \\
4 & -1 \\
\end{pmatrix}
,\ b = \begin{pmatrix}
-5 \\
-9 \\
\end{pmatrix}
</math>
: donc
<source lang="scilab">
-->A = [3 1 ; 4 -1];
-->b = [-5 ; -9];
-->linsolve(A, b)
ans =
2.
- 1.
</source>
: le vecteur
:: <math>\begin{pmatrix}
2 \\
-1 \\
\end{pmatrix}
</math>
: représente la solution ''x'' = 2 et ''y'' = -1.
: Si maintenant on tape
<source lang="scilab">
-->[x0, kerA] = linsolve(A,b)
kerA =
[]
x0 =
2.
- 1.
</source>
: indiquant que la solution est unique.
'''Exemple'''
: S'il y a plus d'inconnues que d'équations :
:: <math>3x + y = 5</math>
<source lang="scilab">
-->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
</source>
: la première commande renvoie une seule solution ; la deuxième commande indique que les solutions sont de la forme
:: <math>\begin{cases}
x = 1,5 - 0,316\,227\,8 \cdot t\\
y = 0,5 + 0,948\,683\,3 \cdot t
\end{cases}</math>
: notons qu'il s'agit d'une solution approchée, puisque la vraie solution est du type :
:: <math>\begin{cases}
x = 1,5 - \frac{1}{3} \cdot t \simeq 1,5 - 0,33\ldots \cdot t\\
y = 0,5 + t
\end{cases}</math>
Si le système n'admet pas de solution, Scilab affiche le message d'erreur : <code>WARNING:Conflicting linear constraints!</code>.
Scilab permet également la résolution symbolique d'un système linéaire avec la fonction <font id="solve1"><code>solve</code></span>
<source lang="scilab">
w = solve(A, c)
</source>
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
: A×''w'' = ''c''.
'''Exemple'''
: Le système d'équation
:: <math>\begin{cases}
x + y = a_1\\
y = a_2
\end{cases}
</math>
: est décrit par
:: <math>\mathrm{A} = \begin{pmatrix}
1 & 1 \\
0 & 1 \\
\end{pmatrix}
,\ c = \begin{pmatrix}
a_1 \\
a_2 \\
\end{pmatrix}
</math>
<source lang="scilab">
-->A = ["1", "1" ; "0", "1"] ; c = ["a1" ; "a2"];
-->solve(A,c)
ans =
!a1-a2 !
! !
!a2 !
</source>
: indiquant que la solution est
:: <math>x = a_1 - a_2</math> et <math>y = a_2</math>.
Scilab propose aussi la [[w:Méthode du gradient conjugué|méthode du gradient conjugué]] avec la fonction <code>conjgrad()</code>. La syntaxe élémentaire est :
<source lang="scilab">
X = conjgrad(A, b)
</source>
=== Système d'équations quelconque ===
Scilab permet la résolution numérique d'un système d'équations quelconques, avec la fonction <code id="fsolve1">fsolve</code>.
Considérons un système de ''n'' équations à ''n'' inconnues (''x''<sub>1</sub>, …, ''x<sub>n</sub>''), de la forme
: <math>\left \{ \begin{matrix}
f_1(x_1, \ldots, x_n) = 0 \\
\vdots \\
f_n(x_1, \ldots, x_n) = 0 \\
\end{matrix} \right .</math>
On définit, en tant que fonction externe (voir ''[[Découvrir Scilab/Programmation#Fonctions|Programmation > Fonctions]]'') la fonction ƒ qui associe à un vecteur '''x'''(''x''<sub>1</sub>, …, ''x<sub>n</sub>'') un vecteur '''v'''(''vx''<sub>1</sub>, …, ''v<sub>n</sub>'') :
: <math>\mathbf{v} = f(\mathbf{x})\text{ : } \left \{ \begin{matrix}
v_1 = f_1(x_1, \ldots, x_n) \\
\vdots \\
v_n = f_n(x_1, \ldots, x_n) \\
\end{matrix} \right .</math>
Le système d'équations consiste donc à résoudre : ƒ('''x''') = '''0'''.
La syntaxe utilisée est :
<source lang="scilab">
[x] = fsolve(x0, f)
</source>
où <code>f</code> est la fonction externe telle que définie ci-dessus et <code>x0</code> 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 <code>f</code> renvoie un vecteur de dimension ''n''.
Scilab utilise un algorithme du type [[w:en:Powell's method|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 :
<source lang="scilab">
[x, v] = fsolve(x0, f)
</source>
Il est recommandé de fournir la jacobienne de ƒ, sous la forme d'une fonction externe ƒ<sub>J</sub> :
: <math>\mathbf{w} = f_\mathrm{J}(\mathbf{x}) = \mathrm{J}\mathbf{x}</math>
où J est la matrice jacobienne :
<source lang="scilab">
[x] = fsolve(x0, f, fJ)
</source>
'''Exemple'''
: En reprenant le système d'équations ci-dessus
:: <math>\begin{cases}
3x+y=5\\
4x-y=9
\end{cases}
</math>
: on écrit donc
<source lang="scilab">
function [y] = f(x)
y = [3, 1 ; 4, -1]*x - [5 ; 9];
endfunction
x0 = [0 ; 0];
solution = fsolve(x0, f)
</source>
: ce qui donne bien (2 ; -1).
Si la fonction est paramétrée, alors
* dans la définition de la fonction, la variable doit venir avant les paramètres ;
* lors de l'appel de <code> fsolve()</code>, la fonction doit être encapsulée avec ses paramètres dans une liste.
'''Exemple'''
: Considérons le système
:: <math>\begin{cases}
a \cdot x + b = 0 \\
c \cdot y + d = 0
\end{cases}</math>
: que l'on veut résoudre pour ''a'' = 1, ''b'' = 2, ''c'' = 3, ''d'' = 4. On écrit donc
<source lang="scilab">
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);
</source>
: ce qui donne comme résultat (-2 ; {{formatnum:-1.333,333,3}}).
=== Équation différentielle ordinaire ===
Une d'équation différentielle ordinaire ''(ordinary differential equation)'', en abrégé EDO ''(ODE)'', peut être résolue de manière numérique avec la fonction <font id="ode1"><code>ode</code></span>. Soit une fonction ''y''(''t'') ; si l'équation différentielle est
: d''y''/d''t'' = ƒ(''t'',''y''),
alors ƒ ayant été définie (fonction externe), la syntaxe pour déterminer ''y''(''t'') est
: <code>''y'' = ode(''y0'',''t0'',''t'',''f'')</code><br />
où
* ''y0'' et ''t0'' sont les valeurs initiales du système, ''y''(''t0'') = ''y0'',
* ''t'' est un vecteur de valeurs pour lesquelles on calcule les solutions, et
* ''y'' est le vecteur de solutions (<font id="plot1"><code>plot(''t'',''y'')</code></span> permet de tracer le graphique des solutions).
La fonction interne <code>ode</code> admet des arguments permettant de résoudre des situations spécifiques :
* <code>"roots"</code> pour rajouter une équation ''g''(''t'',''y'') = 0,
* <code>"discrete"</code> pour calculer de manière récursive ''y''(''k''+1) = ƒ(''k'',''y''(''k'')) à partir d'un état initial ''y''(''k''0).
; Exemple
: On se place sur un segment [0,5;1] ; on veut déterminer les valeurs numériques de la fonction ''y'' vérifiant
:: <math>\frac{\mathrm{d}y}{\mathrm{d}t} = t</math> et <math>y(0) = 0</math>
: on a donc en fait ƒ(''t'',''y'') = ''t'' , t0 = y0 = 0
<source lang="scilab">
-->t=[0.5:0.1:1]; // points de [0,5;1] espacés de 0,1
-->deff("[u] = f(t,y)", "u = t") // définition de la fonction f
-->ode(0,0,t,f)
ans =
0.125 0.18 0.245 0.32 0.405 0.5
</source>
: qui sont bien les valeurs de 1/2×''t''²
; Note
: La fonction ƒ doit être continue et différentiable en tout point. Si la fonction présente des singularité, alors il faut la résoudre sur la première partie D<sup>1</sup> du domaine, puis sur la seconde partie en prenant comme conditions limites les résultats du premier calcul.
; Pour aller plus loin
: Voir {{lien web
| url = http://www.math.univ-metz.fr/~sallet/ODE_Scilab.pdf
| auteur = G. Sallet
| titre = Ordinary Differential Equations with Scilab – WATS Lectures, Provisional notes (Université de Saint-Louis)
| date = 2004 | consulté le = 2018-02-23
| site = Université de Metz
}}
=== Exemple complet ===
[[Fichier:Levage par pivot et appui ponctuel geometrie analytique.svg|vignette|upright=1.7|Levage d'un objet avec un pivot fixe en A et un vérin poussant en B.]]
Considérons un système mécanique simple : une poutre pivote autour d'un point fixe A(''x''<sub>A</sub>, ''y''<sub>A</sub>) ; un vérin fixe pousse sur cette poutre selon un contact ponctuel au point B(''x''<sub>B</sub>, ''y''<sub>B</sub>). La hauteur de la poutre vaut ''e'' = {{unité|200|mm}}, le pivot est implanté à mi-hauteur. Le vérin est implanté à une distance horizontale L = {{unité|1000|mm}}. On a donc ''x''<sub>B</sub> – ''x''<sub>A</sub> = L, et l'altitude ''y''<sub>B</sub> 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 :
: <math>\mathrm{A} \begin{pmatrix} x_\mathrm{A} = 0 \\ y_\mathrm{A} = 0 \end{pmatrix}</math>
: <math>\mathrm{B} \begin{pmatrix} x_\mathrm{B} = \mathrm{L} \\ y_\mathrm{B} \end{pmatrix}</math>
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 :
: (D) : ''y'' + ''ax'' + ''b'' = 0 soit (D) : ''y'' = –''ax'' – ''b''
et l'on veut donc déterminer les paramètres ''a'' et ''b''. Les conditions géométriques sont<ref>voir ''[[w:Propriétés métriques des droites et des plans|Propriétés métriques des droites et des plans]]'' sur Wikipédia.</ref> :
: la droite passe par le point B : ''y''<sub>B</sub> + ''ax''<sub>B</sub> + ''b'' = 0
: la distance de A à (D) vaut ''r'' = ''e''/2 : <math>\frac{y_\mathrm{A} + ax_\mathrm{A} + b}{1 + a^2} - r^2 = 0</math>
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 :
* les inconnues à résoudre sont mises dans un vecteur X : <math>\mathrm{X} \begin{pmatrix} a \\ b \end{pmatrix}</math> ;
* les paramètres du problème sont mis dans un vecteur A : <math>\mathrm{A} \begin{pmatrix} x_\mathrm{A} \\ y_\mathrm{A} \\ x_\mathrm{B} \\ y_\mathrm{B} \\ r \end{pmatrix}</math>.
Le système d'équations devient donc :
: <math>\left \{ \begin{matrix}
\mathrm{A}(4) + \mathrm{A}(3) \times \mathrm{X}(1) + \mathrm{X}(2) = 0 \\
\frac{\mathrm{A}(2) + \mathrm{A}(1) \times \mathrm{X}(1) + \mathrm{X}(2)}{1 + \mathrm{X}(1)^2} - \mathrm{A}(5)^2 = 0
\end{matrix} \right .</math>
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 :
* choisir un point de départ X<sub>0</sub> proche de la solution afin que l'algorithme converge vers le minimum local adapté ; on choisit la droite (A'B), avec
*: A'(''x''<sub>A</sub>, ''y''<sub>A</sub> – ''r'') ;
*: (A'B) : ''y'' = ''y''<sub>A</sub> + (''y''<sub>B</sub> – ''y''<sub>A'</sub>)/(''x''<sub>B</sub> – ''x''<sub>A</sub>)⋅(''x'' – ''x''<sub>A</sub>), donc
*: X<sub>0</sub>(1) = –(''y''<sub>B</sub> – ''y''<sub>A'</sub>)/(''x''<sub>B</sub> – ''x''<sub>A</sub>) = –(A(4) – A(2) + A(5))/(A(3) – A(1)) ;
*: X<sub>0</sub>(2) = (''y''<sub>B</sub> – ''y''<sub>A'</sub>)/(''x''<sub>B</sub> – ''x''<sub>A</sub>)⋅''x''<sub>A</sub> – ''y''<sub>A</sub> = (A(4) – A(2) + A(5))/(A(3) – A(1)) × A(1) – A(2) ;
* vérifier que la solution trouvée est bonne : on remarque que les deux droites possibles sont symétriques par rapport à la droite (AB), il suffit donc de vérifier que la pente de la droite trouvée est supérieure à celle de la droite (AB) :
*: pente<sub>(AB)</sub> = (''y''<sub>B</sub> – ''y''<sub>A</sub>)/(''x''<sub>B</sub> – ''x''<sub>A</sub>) = (A(4) – A(2))/(A(3) – A(1)).
Pour accélérer les calculs, lorsqu'une expression revient plusieurs fois, on la met dans une variable ; on rappelle que seul ''y''<sub>B</sub> = A(4) est variable. On identifie donc :
: facteur<sub>1</sub> = 1/(A(3) – A(1)) ;
: facteur<sub>2</sub> = (–A(2) + A(5)) × facteur<sub>1</sub> ;
: facteur<sub>3</sub> = A(4) × facteur<sub>1</sub> + facteur<sub>2</sub> ;
soit
: X<sub>0</sub>(1) = –facteur<sub>3</sub> ;
: X<sub>0</sub>(2) = facteur<sub>3</sub> × A(1) – A(2) ;
: pente<sub>(AB)</sub> = (A(4) – A(2)) × facteur<sub>1</sub>.
Donc,
* on connaît A(''x''<sub>A</sub>, ''y''<sub>A</sub>, ''x''<sub>B</sub>, ''y''<sub>B</sub>, ''r''){{'}}, la valeur ''y''<sub>B</sub> prenant plusieurs valeurs successives ; on cherche X(''a'', ''b''){{'}} ;
* on calcule les valeurs dérivées facteur<sub>1</sub> = ƒ(A(1), A(3)) et facteur<sub>2</sub> = ƒ(A(2), A(5)) qui sont fixes, et la valeur facteur<sub>3</sub> = ƒ(facteur<sub>1</sub>, facteur<sub>2</sub>, A(4)) qui dépend de ''y''<sub>B</sub> ;
* pour chaque valeur de ''y''<sub>B</sub> (et donc de facteur<sub>3</sub>) :
** on détermine la valeur de départ X<sub>0</sub>,
** on résout le système d'équations, donc on détermine ''a'' et ''b'',
** on vérifie que ''a'' est compatible avec la réalité.
On a donc deux vecteurs de valeurs ''a'' et ''b'', les valeurs ''a''(''i'') et ''b''(''i'') correspondant au déplacement du piston ''y''<sub>B</sub>(''i''), ce qui permet de déterminer l'inclinaison de la poutre pour chaque valeur ''y''<sub>B</sub>(''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.
[[Fichier:Levage par pivot et appui ponctuel isolement poutre.svg|vignette|upright=1.7|Isolement de la poutre, et représentation des caractéristiques connues des forces : on ne connaît que le point d'application de F<sub>A</sub> (flèche brisée), et l'on connaît la direction de F<sub>B</sub> (double flèche avec trait d'axe).]]
Une fois ceci fait, on désire déterminer les effort aux liaisons, c'est-à-dire la force F<sub>A</sub> que l'axe exerce sur la poutre en A, et la force F<sub>B</sub> que le vérin exerce sur la poutre en B. La force F<sub>A</sub> est totalement inconnue :
: <math>\mathrm{F_A} = \begin{pmatrix} x_\mathrm{FA} \\ y_\mathrm{FA} \end{pmatrix}</math>
et on connaît la direction de la force F<sub>B</sub> : c'est le vecteur normal (–''a'', 1). Le vecteur peut donc s'écrire sous la forme :
: <math>\mathrm{F_B} = \mathrm{C}\begin{pmatrix} -a \\ 1 \end{pmatrix}</math>.
Le système a une ''m'', le centre des masses étant en G(''x''<sub>G</sub>, ''y''<sub>G</sub>) ; prenons arbitrairement ''m'' = {{unité|100|kg}}, ''x''<sub>G</sub> = {{unité|700|mm}} et ''y''<sub>G</sub> = {{unité|200|mm}}. L poids P s'écrit p(0, –''mg'') avec ''g'' = {{unité|9.81|m/s<sup>2</sup>}}.
Les équations de la statique s'écrivent
: équilibre des forces : <math>\vec{\mathrm{F}}_\mathrm{A} + \vec{\mathrm{F}}_\mathrm{B} + \vec{\mathrm{P}} = \vec{0}</math>
: équilibre des moments par rapport à A : <math>\mathcal{M}_\mathrm{A}(\vec{\mathrm{F}}_\mathrm{A}) + \mathcal{M}_\mathrm{A}(\vec{\mathrm{F}}_\mathrm{B}) + \mathcal{M}_\mathrm{A}(\vec{\mathrm{P}}) = 0</math>
soit
: équilibre des forces : <math>\left \{ \begin{matrix} x_\mathrm{FA} - \mathrm{C}a + 0 = 0 \\
y_\mathrm{FA} + \mathrm{C} -mg = 0 \end{matrix} \right .</math>
: équilibre des moments par rapport à A :<math> 0 - x_\mathrm{G} \cdot mg + \mathrm{C}\mathrm{L} = 0</math>.
Si l'on définit le vecteur des inconnues
: <math>\mathrm{X_F} = \begin{pmatrix} x_\mathrm{FA} \\ y_\mathrm{FA} \\ \mathrm{C} \end{pmatrix}</math>
et les vecteurs de paramètre
: <math>\mathrm{A_F} = \begin{pmatrix} 1 & 0 & -a \\
0 & 1 & 1 \\
0 & 0 & L \end{pmatrix}</math>
: <math>b_\mathrm{F} = \begin{pmatrix} 0 \\ -mg \\ -x_\mathrm{G} \cdot mg \end{pmatrix}</math>
les équations de la statique s'écrivent sous la forme d'une équation matricielle :
: A<sub>F</sub>X<sub>F</sub> + ''b''<sub>F</sub> = 0.
En résolvant ce système, on peut en déduire les intensités des forces :
: <math>\mathrm{F_A} = \| (x_\mathrm{FA}, y_\mathrm{FA}) \|</math> ;
: <math>\mathrm{F_B} = \mathrm{C} \| (-a, 1) \|</math>.
Concrètement, le programme développé est :
<source lang="scilab">
// 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
</source>
On obtient donc :
<source lang="text">
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.
</source>
== Le problème de la précision et de l'exactitude ==
----
''[[Découvrir Scilab/Calculs élémentaires|Calculs élémentaires]]'' < [[Découvrir Scilab|↑]] > ''[[Découvrir Scilab/OptimisationRésolution d'une fonctionéquations|OptimisationRésolution d'une fonctionéquations]]''
[[Catégorie:Découvrir Scilab (livre)|CalculResolution numériqued'equations]]
| 