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ériques
plot(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(sol,t);           %% 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(sol,t)
       %% y=deval(sol,t)
       %% 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;%% polynomial degre for polyfit
       for k=1:l
              %% first intevals included in [x(1); x(3)]
              x_tmp=x(1:3);
              ind=find((t<x(2)));
              P=polyfit(x_tmp,Y(k,1:3),2);%% polynomial degre=2
              new_y=vec(polyval(P,t(ind))).' ;%% vec and .' for line vectors
              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);%% polynomial degre=3
                     new_y=[new_y vec(polyval(P,t(ind))).']; %% vec and .' for line vectors
              end
              %% first intevals included in [x(1); x(3)]
              x_tmp=x(end-2:end);
              ind=vec(find((t>=x(end-1)))).';
              P=polyfit(x_tmp,Y(k,end-2:end),2);%% polynomial degre=2
              new_y=[new_y vec(polyval(P,t(ind))).'];%% vec and .' for line vectors
              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.