a deval function for octave

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

a deval function for octave

philippe
Hi to all

the syntax to solve a differential equation in matlab has recently
changed and became incompatible with the octave one ! Actually you
should do :

%% solving y''(t)+y(t)=0  y(0)=y'(0)=1 => y(t)=cos(t)+sin(t)
y0=[1;1];T=40;  %%
t=[0:0.1:T];  %%  you choose times t(k)  where solution is evaluated
sol=ode45(@(t,y) [y(2);-y(1)],[0,T],y0);   %%  solving (E)
tt=sol.x;y=sol.y;   %% retrieve tt and y(tt) vectors
clf ;  %% plotting numerical and exact solutions
plot(tt,y(1,:),'xr',t,cos(t)+sin(t),'-b')

but the problem is that you can’t choose times tt computed with ode45 ,
if you want more points you need to do an interpolation between those
computed by ode45 … in matlab this is done by deval function :


y0=[1;1];T=40;%%
t=[0:0.1:T];%%  you choose times t(k)  where solution is evaluated
sol=ode45(@(t,y) [y(2);-y(1)],[0,T],y0);
y=deval(t,sol);
clf ; %% plotting numerical and exact solutions
plot(t,y(1,:),'xr',t,cos(t)+sin(t),'-b')

but deval has not yet been implemented in octave !!! So I’ve tried to
implement my own one (see below), can it be useful to share it here?

sincerely yours ,
Philippe

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=5;
       for k=1:l
              x_tmp=x(1:3);
              ind=find((t<x(2)));
              P=polyfit(x_tmp,Y(k,1:3),K);
              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),K);
              new_y=[new_y polyval(P,t(ind))];
              y(k,:)=new_y;
       end




Reply | Threaded
Open this post in threaded view
|

Re: a deval function for octave

Juan Pablo Carbajal-2
Hi,

I think that the best is to make accessible the interpolation
functions from the ode45 (this deval thingy was done insde the solver
when the user provided a time vector) and use them to implement the
deval functionality.
Could you add an issue to the bug tracker in savahnna?

Thanks
Regards,

Reply | Threaded
Open this post in threaded view
|

Re: a deval function for octave

philippe
Le 25/05/2019 à 00:25, Juan Pablo Carbajal a écrit :
> Hi,
>
> I think that the best is to make accessible the interpolation
> functions from the ode45 (this deval thingy was done insde the solver
> when the user provided a time vector) and use them to implement the
> deval functionality.

I agree that old syntax of ode* looks better, but since matlab has
changed its syntax which now require interpolation with deval "outside"
ode45 you should  had deval to octave for compatibility with matlab.


> Could you add an issue to the bug tracker in savahnna?

I can do it tomorow,  but it's not a bug, it's a feature :-)

best regards,
Philippe


Reply | Threaded
Open this post in threaded view
|

Re: a deval function for octave

philippe
In reply to this post by Juan Pablo Carbajal-2
Le 25/05/2019 à 00:25, Juan Pablo Carbajal a écrit :
> Hi,
>
> I think that the best is to make accessible the interpolation
> functions from the ode45 (this deval thingy was done insde the solver
> when the user provided a time vector) and use them to implement the
> deval functionality.
> Could you add an issue to the bug tracker in savahnna?

https://savannah.gnu.org/bugs/index.php?56393


Philippe



Reply | Threaded
Open this post in threaded view
|

Re: a deval function for octave

Carlo de Falco-2
In reply to this post by philippe
Dear Philippe,

> Il giorno 24 mag 2019, alle ore 11:55, philippe <[hidden email]> ha scritto:
>
> Hi to all
>
> the syntax to solve a differential equation in matlab has recently
> changed and became incompatible with the octave one ! Actually you
> should do :
>
> %% solving y''(t)+y(t)=0  y(0)=y'(0)=1 => y(t)=cos(t)+sin(t)
> y0=[1;1];T=40;  %%
> t=[0:0.1:T];  %%  you choose times t(k)  where solution is evaluated
> sol=ode45(@(t,y) [y(2);-y(1)],[0,T],y0);   %%  solving (E)
> tt=sol.x;y=sol.y;   %% retrieve tt and y(tt) vectors
> clf ;  %% plotting numerical and exact solutions
> plot(tt,y(1,:),'xr',t,cos(t)+sin(t),'-b')
>
> but the problem is that you can’t choose times tt computed with ode45 ,
> if you want more points you need to do an interpolation between those
> computed by ode45 … in matlab this is done by deval function :


This does not seem a recent change to me, AFAIK this syntax has been there for a while (for decades if I am not wrong) ...
If you want your solution to be evaluated a specific points you actually can, just use the syntax (either matlab or octave):

  y0=[1;1];
  T=40;
  t=[0:0.1:T];%%  you choose times t(k)  where solution is evaluated
  [~, y] = ode45(@(t,y) [y(2);-y(1)], t, y0);
  plot (t, y(1,:), 'xr', t, cos(t)+sin(t),'-b')


>
> y0=[1;1];T=40;%%
> t=[0:0.1:T];%%  you choose times t(k)  where solution is evaluated
> sol=ode45(@(t,y) [y(2);-y(1)],[0,T],y0);
> y=deval(t,sol);
> clf ; %% plotting numerical and exact solutions
> plot(t,y(1,:),'xr',t,cos(t)+sin(t),'-b')
>
> but deval has not yet been implemented in octave !!! So I’ve tried to
> implement my own one (see below), can it be useful to share it here?

Your current implementation is not consistent with the numerical method used by ode45 and is
not compatible with matlab, also it does not take into consideration the other ode integrators.

If you are interested in working on a compatible and consistent implementation of deval
that would indeed be a welcom addition to Octave!


> sincerely yours ,
> Philippe

c.