jeudi 18 juillet 2019

une fonction deval pour octave

Résoudre  numériquement une équation différentielle est une fonctionnalité très importante pour un logiciel de calcul numérique comme Scilab, Matlab® ou Octave ... chacun de ces logiciels possède des fonctions pour résoudre ce type de problème, leurs syntaxes sont très proches mais pas forcément compatibles. Si Scilab a choisit de développer sa propre syntaxe  Octave a fait le choix d'être un clone libre de Matlab® et donc de coller au plus près de sa syntaxe. Au contraire Matlab® à tout intérêt  à ne pas être compatible avec un logiciel libre qui proposerait les mêmes fonctionnalités. La question des équations différentielles donne un bon exemple de ce jeu du chat et de la souris entre un modèle propriétaire et son clone libre.






Considérons l'équation différentielle :

$$ (E) ~~~ y''(t)+y(t) =0~~et~~ y(0)=y'(0)=1$$

sa solution est $ y(t)=\cos(t)+\sin(t)$  Pour la calculer numériquement il faut représenter cette équation sous la forme d'un système d'équations différentielles d'ordre 1
$${d\over dt} Y(t)
={d\over dt}\begin{pmatrix}y(t)\\y'(t)\end{pmatrix}
=  \begin{pmatrix}y'(t)\\y''(t)\end{pmatrix}
=  \begin{pmatrix}y'(t)\\-y(t)\end{pmatrix}=f(t,Y(t))
$$
il faut donc  implémenter cette fonction en octave ce qui donne :
 
function y=f(t,x)
   y(1)=x(2);
   y(2)=-x(1); 


on peut alors faire appel à ode45 (ou ode23) pour résoudre numériquement l'équation:

y0=[1;1];T=40;            %%  condition initiale et temps maximal T
t=[0:0.1:T];              %% on choisit les temps t(k) où est calculé y 
[t,y]=ode45(@f,t,y0);     %% résolution numérique de  (E)
clf ;  hold on            %% tracé des solutions exactes et numériques 
plot(t,y(:,1),'xr')
plot(t,cos(t)+sin(t),'-b') 


(on peut remplacer @f par  @(t,y) [y(2);-y(1)]  si on ne veut pas créer de fichier à part pour la fonction f). Hélas cette syntaxe a été changée récemment dans Matlab® et est devenue incompatible avec celle d'Octave! La syntaxe d'octave a donc été modifiée pour accepter la nouvelle syntaxe :


y0=[1;1];T=40;          %%  condition initiale et temps maximal T 
sol=ode45(@f,[0,T],y0); %% résolution numérique de  (E)tt=sol.x;y=sol.y;       %% extraction des résultats
clf;hold on  ;          %% tracé des solutions exactes et numériquesplot(tt,y(1,:),'xr')
plot(t,cos(t)+sin(t),'-b') 

Le problème c'est que maintenant on ne peut plus choisir les points auxquels est évalué la solution, ces points sont choisis par ode45 . Si vous voulez plus de points (pour avoir un tracé précis de la solution, pour calculer une intégrale à partir de la solution, ...) il vous faut faire une interpolation entre les points calculés par ode45 . Dans Matlab® une fonction permet de faire cela de manière très directe, c'est  deval  :


y0=[1;1];T=40;            %% condition initiale et temps maximal T
t=[0:0.1:T];              %% on choisit les temps t(k) où est calculé y 
sol=ode45(@f,[0,T],y0);   %% résolution numérique de  (E)tt=sol.x;y=sol.y;         %% extraction des résultats
y=deval(t,sol);           %% interpolation of solution y aux points de t
clf;hold on  ;            %% tracé des solutions exactes et numériques
plot(t,y(1,:),'xr')
plot(t,cos(t)+sin(t),'-b') 

hélas la fonction  deval n'a pas encore été implémentée dans Octave ce qui le rend incompatible avec la syntaxe de Matlab® !!! Confronté à ce problème j'ai donc codé ma propre fonction deval :

function y=deval(t,sol)
       %% y=deval(t,sol)
       %% sol= solution of ode y'(t)=f(t,y(t))  y(0)=y0
       %%       obtained with sol=ode45(@f,[t0,T],y0)
       %%   t=  vector of times [t(1) ...t(n)]
       %%   y= vector of solution values [y(t(1)) ... y(t(n))]
       %%      interpolated  from sol data
       
       Y=sol.y;
       x=sol.x;
       n=length(t);
       l=size(Y,1);
       y=zeros(l,n);
       K=3;
       for k=1:l
              x_tmp=x(1:3);
              ind=find((t<x(2)));
              P=polyfit(x_tmp,Y(k,1:3),2);
              new_y=polyval(P,t(ind)) ;
              for p=2:length(x)-2
                     x_tmp=x(p-1:p+2);
                     ind=find((t>=x(p))&(t<x(p+1)));
                     P=polyfit(x_tmp,Y(k,p-1:p+2),K);
                     new_y=[new_y polyval(P,t(ind))];
              end
              x_tmp=x(end-2:end);
              ind=find((t>=x(end-1)));
              P=polyfit(x_tmp,Y(k,end-2:end),2);
              new_y=[new_y polyval(P,t(ind))];
              y(k,:)=new_y;
       end



Pour calculer  la solution y(t) en des valeurs t choisies par l'utilisateur le principe est faire une interpolation polynomiale avec la fonction polyfit à partir des points calculés par ode45. La difficulté est de découper la liste de valeurs t en sous-listes t(a:b)  qui correspondent à des intervalles  x(p-1:p+2) calculés par ode45 de telle sorte qu'on puisse faire une interpolation avec des polynômes de degré limité (K=3 ici) pour n'introduire qu'une erreur très faible.

1 commentaire:

Pour écrire des formules mathématiques vous pouvez utiliser la syntaxe latex en mettant vos formules entre des "dollars" $ \$....\$ $ par exemple :
- $\sum_{n=1}^\infty {1\over n^2}={\pi^2\over 6}$ s'obtient avec \sum_{n=1}^\infty {1\over n^2}={\pi^2\over 6}
- $\mathbb R$ s'obtient avec {\mathbb R} et $\mathcal D$ s'obtient avec {\mathcal D}
- pour les crochets $\langle .,. \rangle$ dans les commentaires utilisez \langle .,. \rangle
vous pouvez écrire du html dans les commentaires :
- italique <i> ... </i> gras <b> ... </b>
- lien <a href="http://adresse "> .... </a>