[Octave-bug-tracker] [bug #56393] missing deval function for ode

[philippe Roux]

URL:
  <https://savannah.gnu.org/bugs/?56393>

                 Summary: missing deval function for ode
                 Project: GNU Octave
            Submitted by: rouxph
            Submitted on: Sun 26 May 2019 03:27:19 PM UTC
                Category: Octave Function
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Matlab Compatibility
                  Status: None
             Assigned to: None
         Originator Name: 
        Originator Email: 
             Open/Closed: Open
         Discussion Lock: Any
                 Release: 4.2.2
        Operating System: Any

    _______________________________________________________

Details:

the syntax to solve a differential equation in matlab has recently changed and
became incompatible with the octave one ! Interpolating solution of ode at
specified times is done in matlab by a separate deval function :


y0=[1;1];T=40;
t=[0:0.1:T];%% times t(k)  for solution evaluation
sol=ode45(@(t,y) [y(2);-y(1)],[0,T],y0);  %% solving ode
y=deval(t,sol);  %% interpolating solution
clf ; %% plotting numerical and exact solutions
plot(t,y(1,:),'xr',t,cos(t)+sin(t),'-b')


here is a simple deval function using polyfit/polyval  for interpolation :


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 to this item at:

  <https://savannah.gnu.org/bugs/?56393>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/