Découvrir Scilab/Calcul numérique

Table des matières — Index

---

3. Calcul numérique

---

Vocabulaire

Dans le vocabulaire habituel de la programmation, une procédure est un sous-programme, c'est-à-dire un bout de programme auquel on donne un nom, et que l'on exécute à la demande en « appelant » ce nom. On peut voir cela comme une macro-définition (ou simplement « macro ») : le nom de la procédure est « remplacé » par le contenu du sous-programme. Une fonction est une procédure qui renvoit une valeur.

Dans la documentation officielle de Scilab, les termes function, procedure et macro sont équivalents ; par contre, ils sont utilisés pour les sous-programme définis par l'utilisateur. Les fonctions intégrées dans Scilab, dites built-in, sont appelées « primitives ».

Lorsqu'une primitive a pour argument une fonction, cette fonction est dite « externe » (external). C'est le cas par exemple des primitives de résolution d'équations différentielle ode() ou d'optimisation optim() : on indique à ces primitives quelle est la fonction à utiliser.

Par soucis de simplicité, nous ne feront pas la différence entre les fonctions définies par l'utilisateur et les primitives, sauf lorsque cette nuance est pertinente. Les primitives sont parfois appelées « commandes » ou « instructions » dans le présent document.

Définition de fonctions extérieures

Une fonction extérieure est une fonction numérique définie par l'utilisateur (notez la différence entre « externe » et « extérieure »).

Pour les cas simples, la commande deff permet de définir de nouvelles fonctions pouvant s’exprimer avec les opérateurs déjà prédéfinis (les « fonctions internes » ou « primitives »). On passe deux chaînes de caractères comme paramètre :

• la première indique le nom de la fonction et les variables utilisées en entrée et en sortie, et
• la deuxième indique la formule.

Par exemple, la fonction

ƒ(x) = 2·x

peut se définir par

deff("[y] = f(x)", "y = 2*x")

On peut aussi utiliser l'instruction function() sous la forme

-->function [y] = f(x)
--> …
--> endfunction

la formule, sous la forme y = expression, étant placée à la place des trois points. Avec l'exemple ci-dessus, cela donne :

-->function [y] = f(x)
-->y=2*x
-->endfunction

(voir aussi le chapitre Programmation).

La fonction s'utilise par la suite de la même manière qu'une fonction interne (sur l'exemple ci-dessus) :

-->f(2)
 ans  =

    4.

-->feval(1,f) // calcule f(1)
 ans  =

    2.

Graphe d'une fonction

Il est souvent utile de pouvoir afficher des fonctions, ainsi voici un exemple de la syntaxe Scilab pour afficher la représentation graphique y = ƒ(x) d'une fonction ƒ.

t=[1:0.1:30]; // t est un vecteur des valeurs de 1 à 30 par pas de 0,1

deff("[y] = f(x)", "y = cos(x)./x")  // définition de la fonction f

plot(t, f(t)) // tracé de la fonction

Pour plus de détails voir : Découvrir Scilab/Graphiques et sons#Tracé de fonction.

Gestion des valeurs particulières

La valeur NaN (not a number, résultat d'une opération invalide) peut créer des problèmes pour certaines fonctions (voir Calcus élémentaires > Constantes particulières). La fonction thrownan() permet d'éliminer ces valeurs d'un vecteur ou d'une matrice, et d'indiquer les indice des valeurs non-NaN :

--> x = [0.1 %nan 0.2 0.3]
 x  = 

   0.1   Nan   0.2   0.3

--> thrownan(x)
 ans  =

   0.1   0.2   0.3

--> [sansNaN insensible] = thrownan(x)
 sansNaN  = 

   0.1   0.2   0.3
 insensible  = 

   1.   3.   4.

Gestion du format des nombres

Certaines fonctions ne s'appliquent qu'avec des nombres réels ; et certaines fonctions génèrent des nombres complexes même lorsque le résultat est réel (il s'agit alors d'un nombre complexe dont la composante imaginaire est nulle, du type a + 0⋅i). Il est alors nécessaire de transformer le nombre complexe en nombre réel, avec la fonction real(). Par exemple, nous cherchons les racines du polynôme x² :

--> x = poly(0,"x") // définit le monôme x
 x  = 

  x

--> A = roots(x*x) // calcule les racines de x²
 A  = 

   0. + 0.i
   0. + 0.i

--> isreal(A) // teste si le résultat est réel
 ans  =

  F

--> real(A) // convertit le résultat en nombres réels
 ans  =

   0.
   0.

De la même manière, les nombres entiers sont par défaut représentés par des nombres réels dont la partie décimale est nulle. On peut forcer un format entier avec les commandes intN(), N désignant le nombre de bits utilisé pour représenter le nombre. Par exemple :

--> type(1) // Le nombre « 1 » est de type 1, c'est-à-dire réel ou complexe
 ans  =

   1.

--> type(int64(1)) // Le nombre « int64(1) » est de type 8, c'est-à-dire entier
 ans  =

   8.

Le problème de la précision et de l'exactitude

Scilab est destiné à faire du calcul. La précision, l'exactitude des résultats est donc une préoccupation primordiale.

La première chose est bien évidemment d'utiliser des méthodes, des algorithmes corrects, validés, documentés. Mais la mise en code de ces algorithmes doit également être correcte. Bien évidemment, la méthode choisie doit être retranscrite fidèlement et sans erreur mais il faut aussi prendre en compte certains « effets de bord » dus à la représentation des nombres.

En effet, Scilab calcule par défaut avec des nombre décimaux codés en virgule flottante à double précision selon la norme IEEE754. Cela permet de représenter des nombres dont la valeur absolue est comprise entre 10^-307 et 10³⁰⁸ ; si un calcul donne une valeur inférieure, cela retourne la valeur 0 (erreur de soupassement), et s'il donne une valeur supérieure, cela retourne la valeur Inf ou –Inf (erreur de dépassement). Le nombre de chiffres significatifs est de l'ordre de 16 ; ainsi, si deux nombres diffèrent au 17^e chiffre significatif, ils seront considérés comme égaux (erreur d'arrondi), leur différence sera nulle, leur rapport sera de 1.

Il faut donc s'assurer d'une part que le valeur finale est représentable mais aussi que les résultats des calculs intermédiaires sont eux aussi représentables. Par exemple le calcul de abs(-10^-200) donne bien le résultat 10^-200, alors que sqrt((-10^-200)^2) donne 0, puisque le résultat intermédiaire (-10²⁰⁰)² n'est pas représentable (soupassement). Dans certains cas, un algorithme itératif, cherchant une valeur approchée, peut donner un meilleur résultat qu'un calcul direct utilisant une formule exacte… Par exemple, pour calculer l’hypoténuse d'un triangle rectangle ${\sqrt {a^{2}+b^{2}}}$ (norme euclidienne), on peut utiliser la méthode de Moler et Morrison (1983), ce que fait d'ailleurs la commande Scilab norm().

Le problème de la précision peut parfois se manifester sur des calculs « banals ». Par exemple, le vecteur [-5.1:0.2:5.1] devrait contenir les 51 éléments {–5,1 ; –4,9 ; … ; 4,9 ; 5,1}. Or, les cinquante additions successives de la quantité 0,2, du fait des erreurs d'arrondis, font que le dernier élément calculé est légèrement supérieur à 5,1 (il vaut environ 5,1 + 2⋅10^-15) et donc la liste s'arrête à 4,9. Ce genre de problème peut se gérer avec la commande nearfloat() :

-5.1:0.2:nearfloat("succ", 5.1)

ou bien, pour la génération de vecteurs, avec la commande linspace().

Si l'on doit faire des calculs avec des cas « extrêmes » — nombres très proches, ou bien ayant une valeur absolue élevée ou faible —, on s'intéressera à des algorithmes « robustes ». Ceux-ci font fréquemment intervenir une normalisation des données ou une factorisation pertinente. Dans la mesure du possible, si l'on a le choix entre une fonction « fait maison  » et une fonction Scilab faisant le même travail, on favorisera la fonction Scilab qui est censée être optimisée.

Si des calculs intermédiaires font intervenir des nombres positifs très grands ou très petits, on peut utiliser la fonction logarithme qui a initialement été inventée pour cela. Par exemple, pour calculer

$1-{\frac {365!}{(365-23)!\times 365^{23}}}\simeq 0,507$ avec $365!\simeq 10^{778}\gg 10^{308}$

on peut utiliser la fonction logarithme de la fonction gamma, gammaln() :

1-exp(gammaln(365 + 1)-(gammaln(365 - 23 + 1) + 23*log(365)))

On peut également retravailler la formule en

$1-{\frac {365-22}{365}}\times {\frac {365-21}{365}}\times \ldots \times {\frac {365-0}{365}}$

et utiliser

1 - prod( (343/365):(1/365):(364/365) )

Pour effectuer des tests, Scilab dispose de la constante spéciale %eps qui est égale à la plus petite valeur absolue représentable.

Quoi qu'il en soit, la mise en œuvre d'un algorithme doit toujours s'accompagner de tests avec des valeurs pour lesquelles le résultat est connu. Les résultats de référence peuvent être calculés à la main, obtenus avec un logiciel de référence, ou bien les données peuvent être générées en fonction du résultat attendu.

Par ailleurs, un résultat de qualité fait figurer la précision du calcul.

On pourra se reporter à Baudin 2010a.

Voir aussi

Dans Wikipédia

• Récursivité

Liens externes

Sur Scilab.org

• « Tutorials: Optimization », sur Scilab.org (consulté le 13 février 2019)
• « Tutorials: Statistics », sur Scilab.org (consulté le 13 février 2019)
• [pdf] [Baudin 2010a] [Baudin 2010a] (en) Michaël Baudin, « Scilab is not naive » [PDF], sur Consortium Scilab, 2010 (consulté le 29 décembre 2021)
• [pdf] [Baudin 2010b] [Beaudin 2010b] (en) Michaël Baudin, « Floating point numbers in Scilab », sur forge.scilab.org, novembre 2010 (consulté le 29 décembre 2021)

Autres

• (en) « QL - Quadratic Programming », sur Université de Bayreuth, 2010 (consulté le 24 janvier 2013)
• [pdf] (en) Alan Edelman (MIT), « Open Source and Traditional Technical Computing », sur Scilabtec10, 16 juin 2010 (consulté le 1^er juillet 2017)
• [pdf] anglais David Goldberg, « What every computer scientist should know about floating-point arithmetics », dans ACM Computing Surveys, vol. 23, n^o 1, mars 1991 [texte intégral (page consultée le 2016-08-30)]  ; « verison HTML », sur docs.oracle.com (consulté le 5 janvier 2021)
• [pdf] (en) Willian Kahan, « How Futile are Mindless Assessments of Roundoffin Floating-Point Computation ? », sur Université de Berkeley, 11 janvier 2006 (consulté le 14 mai 2019)

Notes

---

Calculs élémentaires < ↑ > Calcul différentiel et intégral