Divergences entre matlab et octave

Matlab™ ©®  est un des premiers logiciel de calcul numérique créé à la fin des années 80, il reste très utilisé dans le domaine industriel  malgré une syntaxe parfois archaïque et l'existence d'alternatives libres comme GNUoctave, Scilab, ou même Python! Chaque année quand je dois utiliser ce logiciel pour illustrer des cours de calcul numérique, je prépare mes scripts avec GNUoctave et évidement lors des TP je découvre de nouvelles  aberrations "fonctionnalités" toujours plus stupéfiantes  ...

comportement différent entre matlab et GNUoctave
lors de la résolution d'un système linéaire 

une syntaxe archaïque

Ayant commencé par travailler avec Scilab quand j'ai du utiliser matlab j'ai eu le sentiment d'utiliser un langage avec des archaïsmes issus des années 80 :

• la gestion des scripts  en connexion avec le système de fichier est déroutante au début. Un fichier monscript.m   est une commande matlab monscript  (youpi c'est facile) mais seulement s'il est dans le répertoire courant (pourquoi ?!?)
•  impossible d'aller exécuter un script d'un autre répertoire sauf à donner un statut particulier à se répertoire comme ceux contenant l'installation de matlab
• pareil pour les définitions de fonctions  qui doivent être dans le répertoire courant avec le même nom que la fonction (plus l'extension .m)  donc une fonction "foo"  doit obligatoirement être dans le fichier foo.m et évidement on ne peut pas (facilement)  charger dans l'espace de travail des fonctions venant d'un autre répertoire  ...

Autre exemple  avec les variables protégées. En Scilab on a accès aux constantes fondamentales $e,\pi,i$  sous forme de variables  %e, %pi ,%i  "protégées"  (c'est la signification du '%')  c'est à dire qui ne peuvent être modifiées par l'utilisateur (ça semble logique et prudent). Hé bien en matlab ces variables e,pi,i  ne sont pas protégées !?! Ca veut dire que si vous écrivez pi=3  quelque part vous aurez perdu la valeur numérique de pi utilisée à tout plein d'endroit dans vos scripts ! Bon on se dit que c'est pas grave on se souviendra pour pi, mais ça peut être plus pernicieux, par exemple si vous avez une boucle comme celle ci-dessous :

for i=1:10
            disp(2^i)
end
exp(i*pi/4)

hé bien exp(i*pi/4) ne donnera pas 0.7071 + 0.7071i  par ce que "i" à été remplacé par 10 à la fin de la boucle 😤

Autre archaïsme : la gestion des polynômes. Un polynôme $P(X)=a_n X^n+...+a_1 X +a_0$  est naturellement représenté par le vecteur le ses coefficients [a_n ... a_1 a_0], un choix simple apriori (à condition de se rappeler de l'ordre des coefficients) mais pas du tout pratique quand on y réfléchi : 

• comment faire le produit de deux polynômes p et q ? pas en faisant p*q ni même p.*q non il faut une fonction spéciale conv(p,q) (hé oui il y a une convolution cachée derrière ce produit)
• comment faire la somme de deux polynômes ? Facile on fait p+q, mais il faut qu'ils aient le même degré sinon ça marche pas 😫
• comment afficher l'expression du polynôme appliqué à une variable 'x' ? Facile il y a une fonction pour ça poly2sm(p,'x') mais elle est dans le package 'symbolic' qui est payant ! Ou sinon il y a la commande polyout(p,'x') mais ça c'est uniquement avec GNUoctave évidement ...

>> p=poly([1 2]) % p= coefficients de (x-1)*(x-2)
p =
   1  -3   2
>>polyout(p,'x') % affichage en octave seulement !
1*x^2 - 3*x^1 + 2
>>conv(poly(1),poly(2)) % produit  (x-1)*(x-2)
ans =
   1  -3   2

comparé à Scilab ça me semble moins ergonomique. Avec cet autre logiciel  on a un type "polynôme" à part qui permet de les additionner/multiplier  avec les opérateur +/* tout en affichant le résultat sous une forme lisible directement:

-->p=poly([1,2],'x')  // p=(x-1)*(x-2) en scilab
 p  =
  2 -3x +x^2
-->x=poly(0,'x')  // monôme x
 x  =
  x

-->(x-1)*(x-2)  // produit en scilab
 ans  =
  2 -3x +x^2

Toutes ces complication ne sont pas si graves mais au fil des années il s'y ajoute des choses que je trouve plus graves ...

Résolution des systèmes linéaires

Depuis quelques années matlab  procède à des modifications de syntaxe parfois incompréhensibles comme si on voulait casser la compatibilité avec GNUoctave  qui deviendrait un vrai concurrent 🤔. Par exemple  la constante e=   2.718281828459046  a disparu du jour au lendemain de matlab, on peut toujours récupérer cette valeur avec la commande exp(1) mais cela peut vous obliger à modifier vos anciens scripts si vous voulez qu'ils fonctionnent avec matlab. 

Plus grave quand un changement modifie en profondeur le comportement d'une commande aussi fondamentale que celle chargée de résoudre un système linéaire d'équation. Cette résolution est très pratique en matlab: il suffit de créer la matrice M et le vecteur y  représentant matriciellement le système M*x=y d'inconnue x et de taper x=M\y ou x=linsolve(M,y). C'est simple et ça applique la méthode de Gauss  au système pour renvoyer une solution dans tous les cas :

• l'unique solution si c'est un système de Cramer (autant d'équations indépendantes que d'inconnues)
• une des solutions si c'est un système sous-déterminé ( moins d'équations indépendantes que d'inconnues)
• une solution au sens des moindres carrés si c'est un système sur-déterminé (plus d'équations indépendantes que d'inconnues ou équations incompatibles) 

ces commandes sont faciles à manipuler, par exemple pour un système de Cramer :

>> M1=[5   7   1;   5   7   3;   5   1   5]%% obtenu avec  randi(9,[3,3])
M1 =

   5   7   1
   5   7   3
   5   1   5

>> det(M1)%% système de Cramer car det(M1) non nul, sinon recommencer
ans = 60
>> x1=[ 7; 4 ; 2]%% obtenu randi(9,[3,1])
x1 =

   7
   4
   2

>> y1=M1*x1 %%
y1 =

   65
   69
   49

>> M1\y1 %% on retrouve bien  x1
ans =

   7
   4
   2

>> linsolve(M1,y1) %% même résultat
ans =

   7
   4
   2

De même est enlevant une équation on obtient un système sous déterminé :

>> M2=M1; M2(2,:)=[], %% enlever équation n°2
M2 =

   5   7   1
   5   1   5

>> y2=M2*x1 %% recalculer y
y2 =

   65
   49

>> x2=M2\y2 %% nouvelle solution  différente de x1
x2 =

   5.643322475570033
   4.798045602605867
   3.197068403908795

>> erreur=norm(M2*x2-y2)%%  vaut bien 0 (aux erreurs numériques près)
erreur = 2.842170943040401e-14

mais si on avait modifié l'équation n° 2 pour la rendre insoluble on aurait une solution au sens des moindre carrés :

>> M3=M1; M3(2,:)=0,y3=y1 %% équation n° 2 du type 0=4
M3 =

   5   7   1
   0   0   0
   5   1   5

y3 =

   65
   69
   49
>> x3=M3\y3 %   solution  au sens des moindres carrés
warning: matrix singular to machine precision
warning: called from
    /tmp/octave_nYpLlc.m at line 1 column 3

x3 =

   5.643322475570034
   4.798045602605862
   3.197068403908794

>> erreur=norm(M3*x3-y3) %% non nul mais minimal
erreur = 69

En plus on a un warning signalant que la solution n'en est pas exactement une! Hélas ce comportement a été cassé dans les versions récentes de matlab sauf à utiliser un paramètre optionnel, qui doit être écrit sous une forme absolument  incompréhensible pour le débutant (option.RECT=true)  🙀 :

>> x3=M3\y3 %   ne trouve pas la solution
Warning: Matrix is singular to working precision.

x3 =

  -Inf
   Inf
   Inf

>> linsolve(M3,y3) % ne trouve pas la solution
Warning: Matrix is singular to working precision.

ans =

  -Inf
   Inf
   Inf

>> option.RECT=true , linsolve(M3,y3,option) % trouve la solution !!!

option =

  struct with fields:

    RECT: 1

Warning: Rank deficient, rank = 2, tol =  4.710277e-15.

ans =

   9.266666666666667
   2.666666666666662
                   0

donc c'est quand même possible avec matlab mais ça complique énormément les choses.

Précision des calculs 

Arrivé là je pensais avoir tout vu, mais non, il était possible de faire pire : affaiblir la précision des calculs ! Je m'en suis rendu compte sur une des choses que je trouve la mieux faite dans matlab : le calcul numérique d'intégrales doubles. Après avoir longtemps multiplié les fonctions d'intégration numériques avec des noms en quad* aussi obscurs que possibles, le code s'était enfin clarifié avec trois fonctions :

• integral(f,a,b) pour calculer des intégrales simples  $\int_a^b f(x) dx$
• integral2(f,a,b,c,d) pour calculer des intégrales doubles $\int_a^b \int_c^d f(x) dy dx$
• et évidement integral3 pour des intégrales triples
  

Ce qui est vraiment bien pensé c'est que $c$ et $d$ peuvent être des fonctions de $x$ pour traiter le cas d'un domaine non-rectangulaire (c'est rare que je dise du bien de matlab!), par exemple l'intégrale :

$$I=\iint_\Omega {\ln(2+x)\over \arctan(x-y)}dx dy ~~sur~~ \Omega=\{(x,y)\vert 0\leq x\leq 2~et~ -x\leq y\leq\sin(x)\}$$

se calcule avec les commandes :

>>f=@(x,y) log(2+x)./atan(x-y);
c=@(x) -x;
d=@(x) sin(x) ;
integral2(f,0,2,c,d)

ans = 7.453230816525963

 Mais la question quand on fait un tel calcul est de savoir a priori quelle est sa précision (alors qu'on ne connaît pas la valeur exacte de l'intégrale). C'est le genre de choses qui doit être écrit dans la documentation (faire doc integral2)  et on trouve ici un paragraphe qui explique que la précision par défaut est de 1e-10 :

AbsTol Define the absolute error tolerance for the quadrature. The default value is 1e-10 (1e-5 for single).

ce paragraphe est tiré de octave mais on a la même spécification pour matlab dans la documentation, pourtant quand on fait le calcul de $I$  on obtient:

• 7.453231792656116  avec matlab
• 7.453230816525963 avec  octave

les deux résultats ne coïncide qu'à 1e-6 près, mais alors lequel a une précision de 1e-10 ???  Testons en imposant une précision plus grande  avec la commande integral2(f,0,2,c,d,'AbsTol',1e-12)  et on obtient :

• 7.453230816525963 avec matlab
• 7.453230816525963 avec octave

donc c'est bien octave qui donne par défaut une précision de 1e-10 matlab n'atteignant que 1e-6 😭

J'avais déjà noté des différences de résultats dans la résolution numériques d'équations différentielles entre matlab/GNUoctave avec l'équation 

$$ \begin{array}{rcl}x'(t)&=&-4y(t)\\y'(t)&=&\sin(x(t))\\x(0)&=&0\\y(0)&=&-1\end{array}$$

qui se résout dans les langage avec la commande ode45 (méthode de Runge-Kuta d'ordre 4) :

t=[0:0.01:15];u0=[0;-1]; % données du problèmes
f=@(t,u) [-4*u(2); sin(u(1))]; % modélisation de l'équation
[t,u]=ode45(f,t,u0); % résolution
x=u(:,1);y=u(:,2); % extraction des solutions
%%tracé des solutions
clf
subplot(1,2,1)
hold on, grid on
plot(t,x,'-b','linewidth',2)
plot(t,y,'-r','linewidth',2)
legend('x(t)','y(t)')
title("x'(t)=-4y(t), y'(t)=sin(x(t)), x(0)=0, y(0)=-1")
subplot(1,2,2)
hold on,grid on
plot(x,y,'-g')
title('portait de phase (x(t),y(t))')

quand on compare les résultats entre les deux logiciels on obtient les figures suivantes :

divergence entre les solutions d'une équation différentielle pour matlab et GNUoctave

les solutions du système  sont forcément périodiques et se comportement apparaît bien sur GNUoctave, mais on voit que la simulation sur matlab (à droite)  conduit à un changement de comportement près du point de coordonnée $x=\pi,y=0$ (point singulier). C'est en comparant les solutions d'étudiants (travaillant avec matlab) aux miennes (obtenues avec GNUoctave) que j'ai découvert ces exemples, il doit sûrement y en avoir beaucoup d'autres encore plus préoccupant dans un cadre industriel.

Mais il y a encore plus surprenant : si on remplace ode45 par ode23, qui implémente la méthode de Runge Kutta d'ordre 2 (donc moins précise a priori que ode45) matlab donne alors une solution périodique comme le fait encore GNUoctave ?!? A noter que dans scilab on a une fonction ode qui choisit elle-même l'algorithme de résolution, sauf si l'utilisateur indique le choix de méthode par un paramètre optionnel. C'est certainement une approche plus ergonomique que d'imposer à l’utilisateur de choisir entre ode23 et ode45 (qui ne sait pas forcément qu'elle est le meilleur choix par défaut).

Conclusion

Aujourd'hui matlab est  moins stable, moins ergonomique et moins précis que GNUoctave  alors qu'il s'agit d'un logiciel commercial, pourtant il reste en position dominante !?! Comment est-ce possible ? Essentiellement grâce à l'offre de toolboxes (boîtes à outils) reposant sur matlab et qui permettent souvent de contrôler des matériels avec des protocoles propriétaires. L'offre de toolboxes de GNUoctave a beaucoup progressée ces dernières années mais cela doit rester un frein pour des entités qui utilisent matlab  uniquement par rapport à une de ces toolboxes. Le seul élément qui tend à changer la donne c'est le prix exorbitant des licences matlab dans un contexte économique contraint ...