bugGNU Octave - Bugs: bug #56393, missing deval function for ode

 
 

bug #56393: missing deval function for ode

Submitter:  philippe Roux <rouxph>
Submitted:  Sun 26 May 2019 03:27:19 PM UTC
   
 
Category:  Octave Function Severity:  1 - Wish
Priority:  3 - Low Item Group:  Feature Request
Status:  Confirmed Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Tue 28 May 2019 06:19:24 AM UTC, comment #8: 

rik, to add a bit more detail, according to the ode documentation, one of the mandatory fields in the "sol"
structure is the name of the solver used to compute the solution.
The other mandatory fields are only the integration points and solution values "x" and "y".

In addition to the mandatory fields, extra data is stored in "sol". What data is actually there does not seem to be documented and is probably dependent on the solver used.

by running a test with ode45 on matlab 2014a, some fields that seem relevant for the use by deval are

+varargin+
sol.extdata.{odefun,options,varargin}
-varargin-

which appear to be the original inputs to the call to ode45
and

+varargin+
sol.idata.f3d
-varargin-

which looks like the derivatives in each integration interval

all this is not documented so I don't think we should rely on the existence of any of these fields nor copy the exact field names.


Carlo de Falco <cdf>
Group Member
Mon 27 May 2019 10:54:26 PM UTC, comment #7: 

I have checked a simple ODE with Matlab and interp1 does not reproduce the output of deval (with various methods such as "linear", "cubic", and "pchip").  Hence, I think Carlo's ideas are correct: the 'sol' structure must contain the interpolator from the original ODE function and that is used to produce the requested output.

Rik <rik5>
Group administrator
Mon 27 May 2019 09:55:14 PM UTC, comment #6: 

Re: comment #1

I can check with matlab versions back to 2014a (released Feb 11 2014) and there seems to have been no syntax change since then, so I don't see what is meant here by "the syntax to solve a differential equation in matlab has recently changed"

Carlo de Falco <cdf>
Group Member
Mon 27 May 2019 09:44:58 PM UTC, comment #5: 

Re: comment #3

The latest version of matlab I have access to is 2018b and the code I posted works there.

Anyway that behaviour is still documented in the current web documentation for the latest Matlab version so I assume it is still supported.

Re: comment #4

IIRC there was someone looking into how this is handled in Matlab at the time when we ported ode45 and ode23 from odepkg into Octave.

AFAIR the basic idea is that


[t,y] = ode45(odefun,t,y0)


and


sol = ode45(odefun,t,y0)
y = deval(sol,t)


should return exactly the same results.
for this to work Matlab stores into the "sol" structure
all the information needed to use the same interpolator
as is used during the integration.

In particular for ode45 and ode23 the solution and derivatives
at the intermediate stages of RK5 and RK3 method are stored and
the appropriate Hermite interpolation methods are used.

Carlo de Falco <cdf>
Group Member
Mon 27 May 2019 05:51:59 PM UTC, comment #4: 

The Matlab documentation speaks of the output being interpolated, but does not seem to describe how smooth the result should be.  Is it linear, cubic?  The second output of the deval function is a list of first derivatives at the desired points x which suggests that the output is only linear.

In that case, rather than a for loop with polyfit/polyval (which will be slow because of interpreted loops), it might be better to use interp1 and create a 'PP' structure which can then be evaluated in a single operation.

Another possibility would be to use griddatan which relies on a Delaunay triangulation and might be faster, although not if you need the first derivatives as a potential output.

Rik <rik5>
Group administrator
Mon 27 May 2019 02:21:28 PM UTC, comment #3: 

Ok, great.I thought that code wouldn't work in matlab anymore.

Juan Pablo Carbajal <juanpi>
Group Member
Mon 27 May 2019 09:31:07 AM UTC, comment #2: 

@JPi I don't see any incompatibility just the function deval is missing I think the group is correct.

BTW as pointed out on the mailing list, the use of deval can
be avoided with the following code:


 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')


that works with both matlab and octave

Carlo de Falco <cdf>
Group Member
Mon 27 May 2019 07:36:33 AM UTC, comment #1: 

Hi mike,

Why did you change the group out of "Matlab compatibility"?
Matlab code using "deval" will fail in GNU Octave making it ML incompatible.

Regards,

Juan Pablo Carbajal <juanpi>
Group Member
Sun 26 May 2019 03:27:19 PM UTC, original submission:  

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




philippe Roux <rouxph>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by cdf (Posted a comment)
  • -email is unavailable- added by juanpi (Posted a comment)
  • -email is unavailable- added by rouxph (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 7 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-05-27 rik5 StatusNeed Info Confirmed
    2019-05-27 rik5 StatusNone Need Info
    2019-05-26 mtmiller Priority4 3 - Low
    2019-05-26 mtmiller Severity3 - Normal 1 - Wish
        Priority5 - Normal 4
        Item GroupMatlab Compatibility Feature Request
        Release4.2.2 dev

    Back to the top

    Powered by Savane 3.13-3230.
    Corresponding source code