XHTML
Le langage HTML/Introduction/Exercice 1
Découvrir Scilab/Calcul numérique
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 codice_1 ou d'optimisation codice_2 : 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 codice_3 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 :
Par exemple, la fonction
peut se définir par
deff("[y] = f(x)", "y = 2*x")
On peut aussi utiliser l'instruction codice_4 sous la forme
-->function [y] = f(x)
--> endfunction
la formule, sous la forme codice_5, é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
Gestion des valeurs particulières.
La valeur NaN ("", 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 codice_6 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 codice_7. 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 codice_8, "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 10308 ; 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 17e 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 codice_9 donne bien le résultat 10-200, alors que codice_10 donne 0, puisque le résultat intermédiaire (-10200)2 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 formula_1 (norme euclidienne), on peut utiliser la méthode de Moler et Morrison (1983), ce que fait d'ailleurs la commande Scilab codice_11.
Le problème de la précision peut parfois se manifester sur des calculs « banals ». Par exemple, le vecteur codice_12 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 codice_13 :
-5.1:0.2:nearfloat("succ", 5.1)
ou bien, pour la génération de vecteurs, avec la commande codice_14.
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
on peut utiliser la fonction logarithme de la fonction gamma, codice_15 :
1-exp(gammaln(365 + 1)-(gammaln(365 - 23 + 1) + 23*log(365)))
On peut également retravailler la formule en
et utiliser
1 - prod( (343/365):(1/365):(364/365) )
Pour effectuer des tests, Scilab dispose de la constante spéciale codice_16 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 à .
Notes.
"Calculs élémentaires" < ↑ > "Calcul différentiel et intégral"