bugGNU Octave - Bugs: bug #63162, Inaccurate "Events"...

 
 

bug #63162: Inaccurate "Events" location in ode45, ode23(s)

Submitter:  Ken Marek <m5318>
Submitted:  Wed 05 Oct 2022 05:07:44 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Mon 21 Nov 2022 04:51:09 PM UTC, comment #29: 

Yes, I think so. tests pass and doesn't seem to impact anything outside of these functions. can always reopen if something pops up.

Nicholas Jankowski <nrjank>
Group Member
Mon 21 Nov 2022 04:02:05 PM UTC, comment #28: 

OK to close this report?

John W. Eaton <jwe>
Group administrator
Wed 16 Nov 2022 09:37:14 PM UTC, comment #27: 

Re: comment #24
The answers to your questions are both "yes": STEPPER has always been called before runge_kutta_interpolate, and STEPPER always returns a set of valid k_vals. runge_kutta_interpolate is called either directly from integrate_adaptive (which calls STEPPER and therefore has k_vals), or, with this bug fix, from ode_event_handler (which is called from integrate_adaptive, and therefore has k_vals).

Ken Marek <m5318>
Fri 11 Nov 2022 03:51:31 PM UTC, comment #26: 

Sorry for that and apologies for the bad push. Misread the patch split as referring to being across the separate ode bug reports.

Nicholas Jankowski <nrjank>
Group Member
Fri 11 Nov 2022 02:56:13 PM UTC, comment #25: 
John W. Eaton <jwe>
Group administrator
Fri 11 Nov 2022 02:07:07 PM UTC, comment #24: 

I backed out changeset 88fff8521d7.

I asked for the changes to be made separately, not just for the purpose of reviewing them here, but also because it makes sense to me to have changes with smaller scope for people who look at the history later.

In rkinterp_bug63163_20221006.patch, I see that the DT argument to runge_kutta_interpolate is never used.  I also see that if we know that k_vals never empty, then we don't need to call FCN, so it and ARGS become unused arguments as well.  Is it guaranteed that k_vals will never be empty because it is obtained by calling STEPPER in integrate_adaptive prior to calling runge_kutta_interpolate?  And STEPPER has has already evaluated the user-supplied function so there should be no need for runge_kutta_interpolate to also call it?

John W. Eaton <jwe>
Group administrator
Fri 11 Nov 2022 03:40:52 AM UTC, comment #23: 

(apologies, when combining the patches I forgot to put your name on the commit)

Nicholas Jankowski <nrjank>
Group Member
Fri 11 Nov 2022 02:58:56 AM UTC, comment #22: 

ok. i combined the patches, made a few minor text tweaks, and pushed to default as

https://hg.savannah.gnu.org/hgweb/octave/rev/88fff8521d76

marking as Ready for test

Nicholas Jankowski <nrjank>
Group Member
Thu 10 Nov 2022 11:28:39 PM UTC, comment #21: 

Yes, this one is ready to go. Please apply Comment #16 patch, then Comment #17 patch.

Based on the feedback I received here, I have not made any changes to how extra function parameters are handled, except to the extent required to fix the Event finding. I think that particular issue is low priority, as at this point it is mostly a code cleanup project, rather than functionality (unless the feature is deprecated at some point). If I decide to work on it at a later time, I will open a new bug report.

Ken Marek <m5318>
Thu 10 Nov 2022 11:20:11 PM UTC, comment #20: 

ken, sorry for leaving this for a month. was this patch ready to be applied and tested? or was there anything else to do for this item?

Nicholas Jankowski <nrjank>
Group Member
Mon 10 Oct 2022 05:56:15 PM UTC, comment #19: 

That should have read __ode15__.cc in comment #18. One of the underscores got eaten up as a italic format marker, so I'm reposting the corrected filename.

Ken Marek <m5318>
Mon 10 Oct 2022 05:47:51 PM UTC, comment #18: 

It looks like the y cell input for ode_event_handler that I took out was once used by ode15i and/or ode15s. These both currently call the compiled function ode15 (I think from _ode15_.cc), and not the files that I have been amending. I'm not sure if they are also responsible for the unused items in runge_kutta_interpolate; I mostly checked that they do not call the function.

Ken Marek <m5318>
Thu 06 Oct 2022 10:34:14 PM UTC, comment #17: 

The extra function inputs turned out to be a bigger beast than I originally thought, so I have gone back and just implemented the improved Event finding, while leaving the extra input feature intact (there was still a bit of cleanup to make it work with the bug fix). Patch attached. Please apply after the patch in comment #16, in which I accidentally mistyped the bug number in the file name.

Dealing with the extra function inputs can now be done in a different bug report if desired. I did add a FIXME note suggesting a warning for the undocumented feature.

(file #53811)

Ken Marek <m5318>
Thu 06 Oct 2022 08:21:48 PM UTC, comment #16: 

I started making the changes in reverse order to make them separable. Here is the the one that cleans up unused/unnecessary stuff from runge_kutta_interpolate. As far as I have tested, there is no resulting change in behavior for any of the odexx functions.

(file #53810)

Ken Marek <m5318>
Thu 06 Oct 2022 10:35:38 AM UTC, comment #15: 

@m5318, sorry for comment #14 I had tried to apply the patch before pulling changeset 31263:449ed6f427cb with that in place
it works.

Carlo de Falco <cdf>
Group Member
Thu 06 Oct 2022 10:13:35 AM UTC, comment #14: 

@5318, how are you generating your patches?
They never worked for me on the default branch ...

Carlo de Falco <cdf>
Group Member
Thu 06 Oct 2022 09:54:30 AM UTC, comment #13: 


comment #8:

> It depends what you mean by leaving it as is. To do the Events root finding I need to send a handle of the complete Events function to the root solver, so if there are extra inputs I still have to wrap the Events function in some kind of anonymous function or something in ode_event_handler.m



If I understand correctly the use of functions with extra parameters is no longer documented in Matlab.

The Matlab docs were suggesting the use of an anonymous function do bind the extra parameters in the calling code prior to invoking
the ode solvers.

Do we still need to support this feature in Octave even if Matlab has deprecated it?

Carlo de Falco <cdf>
Group Member
Thu 06 Oct 2022 01:13:17 AM UTC, comment #12: 

Posting them all here is fine.

John W. Eaton <jwe>
Group administrator
Wed 05 Oct 2022 10:39:38 PM UTC, comment #11: 

I'll try to look at this again tomorrow. Maybe it won't be as much of a hassle to separate as I think.  If I do three patches, do you want separate bug reports for all three issues, or all of them here?

Ken Marek <m5318>
Wed 05 Oct 2022 09:53:02 PM UTC, comment #10: 

Option 1 would let me separate things out into three separable problems: 1) Events accuracy; 2) Extra function parameters; and 3) Runge-Kutta interpolation function. They could be three separate patches.

However, each patch would have to modify all the same function definitions and function calls in all the same files, with some calls happening several times per file, so it would be maybe 3x the work I have already done, plus (I think) waiting for each patch to push before the next one goes in, since they all modify the same lines. Implementing my comment #6 idea (Option 3) would just mean keeping all my changes so far, plus adding a bit of the same code to each of ode23/23s/45, which is a lot easier at this point (assuming everything is as easy as I think, which is never guaranteed). That is why it is my preference.

To your other comments:
Yes, I think your feval improvements are probably good. I was just trying to add clarity with the less efficient code, and I wasn't sure that self-referring to "fcn" in the assignment and feval wouldn't cause a problem, so I played it safe for the moment.

If I do all the function wrapping in the base odexx file, I'll want to store funarguments so I don't have to do all the input testing again for Events and OutputFcn, where they are also used.

Ken Marek <m5318>
Wed 05 Oct 2022 09:28:26 PM UTC, comment #9: 

It would be good to have separate patches for each problem you're trying to fix.  If that's not possible, or would mean a lot of extra work at this point, then we can work with one larger patch.  But in the future, it would be best to try to limit the scope of changes and not lump too many fixes in one change.

The change you proposed in comment #6 is OK, but I think you could just write


fcn = @(t, y) feval (fcn, t, y, varargin{2:end});


If you do need to have "funarguments" defined some other reason, then you need to use


fcn = @(t, y) feval (fcn, t, y, funarguments{:})


to pass the extra arguments individually to the original user-supplied function instead of collected in a single cell array argument.

Also, I assume "feval" is needed because the original function could be passed as a function handle, inline function, or character array that names the function so that the simple function call syntax "fcn (t, y, ...)" won't work.

John W. Eaton <jwe>
Group administrator
Wed 05 Oct 2022 09:00:02 PM UTC, comment #8: 

It depends what you mean by leaving it as is. To do the Events root finding I need to send a handle of the complete Events function to the root solver, so if there are extra inputs I still have to wrap the Events function in some kind of anonymous function or something in ode_event_handler.m (Option 1), or else copy and modify the entirety of fzero.m for use as a helper function in ode_event_handler.m (Option 2, a real pain but possibly makes some minor efficiency improvements). I can do one of these if you want.

I think it would actually be easier to implement my suggestion from comment #6 to maintain current functionality (Option 3), than to make everything work with as much of the old code as possible, in which case I think I'd still have code modifications in nearly all of the same places.

I can make Option 1 or 3 happen reasonably quickly, if you let me know your preference. Option 2 will take a while.

There's also the other changes I mentioned:

  • Cleaning up runge_kutta_interpolate.m, which I can mostly put back as before or leave fixed; and
  • Removing the y cell input from ode_event_handler.m. This one has to go, as there's no way to accommodate the feature without knowing the underlying model of the non-existent file that would theoretically use it.

So I guess let me know what you want for runge_kutta_interpolate.m as well.

Ken Marek <m5318>
Wed 05 Oct 2022 08:22:15 PM UTC, comment #7: 

For this bug report, can you create a patch that changes only what is necessary to fix the incorrect handling of the Events function and leave the extra function parameters as they are now?

John W. Eaton <jwe>
Group administrator
Wed 05 Oct 2022 07:57:31 PM UTC, comment #6: 

Ok, I think I can make that happen. Just to clarify, this is what you want:

  • Additional parameters passed to ode45, etc. will continue to be passed to fcn, as well as to OutputFcn and Events if they exist, but we will now generate a warning that this functionality is deprecated.


Rather than passing along all the arguments explicitly, couldn't we just redefine the functions in the first script? I.e., in ode45, something like:

temp_fcn = fcn;
funarguments = {varargin{2:numel (varargin)}};
fcn = @(t,y) feval(temp_fcn,t,y,funarguments);

and so on for OutputFcn and Events? Then everything that was given those function handles afterward wouldn't need to worry about whether or not there are extra input parameters.

Ken Marek <m5318>
Wed 05 Oct 2022 07:45:44 PM UTC, comment #5: 

It's probably not necessary to document the undocumented Matlab feature, but I do think it would be good to warn that it will actually be removed in a future version.

John W. Eaton <jwe>
Group administrator
Wed 05 Oct 2022 07:41:47 PM UTC, comment #4: 

i'm looking through the bug #63063 patch now. Making some cosmetic edits at the moment, will likely push it soon if all looks good. So go ahead and work on this one assuming that will be implemented in a form more or less like you have it.  I added a dependency to the report.

Nicholas Jankowski <nrjank>
Group Member
Wed 05 Oct 2022 07:38:33 PM UTC, comment #3: 

Thanks for the information. I'm happy to back out those changes and keep sending along extra arguments to everything if that's what you want; in that case the documentation should also be updated to let users know about the feature, even if it's being deprecated.

Patch file is attached. You will need to first apply the patch from Bug #63063, Comment #21.

I also forgot to mention some changes to runge_kutta_interpolate.m. I could not find any function calling it that did not send a valid k_value matrix, so the function evaluation backups were unneeded. I removed them and the corresponding function input arguments. Otherwise I think all the changes not already mentioned are just documentation.

Since I'm posting a patch, let's do the results too:

octave:130> test_file
Elapsed time is 0.089494 seconds. % cleared memory
Elapsed time is 0.0599439 seconds. % subsequent run
y2 function errors:
t error of event: 6.6613e-16
y error of "true" root: 4.4409e-16
y error of Event root, as reported by event: 2.0123e-15
y error of Event root, from y2(t_event): 3.1086e-15
y7 function errors:
t error of event: 1.8769e-08   5.246e-07
y error of "true" root: 3.1086e-15  4.2633e-14
y error of Event root, as reported by event: 2.5275e-15  4.2077e-14
y error of Event root, from y7(t_event): 9.9752e-08  1.8531e-05

You can see that the Octave run time increases a bit (19% first run, or 5% subsequent for this test case), but we also get results within an order of magnitude of the Matlab results, which is a huge improvement in accuracy. The patch also includes updated tolerances on the file tests for ode23/45 files - Events tests. The test cases now improve accuracy by about 2 decimal places. Ode23s didn't have those tests, so I added them as well.

(file #53797)

Ken Marek <m5318>
Wed 05 Oct 2022 06:45:09 PM UTC, comment #2: 

Matlab did not always have anonymous functions to capture additional parameters and I'm fairly certain that older versions of ode45 in Matlab accepted extra parameters to pass to user-defined functions.  If that feature is really no longer present in these functions in Matlab, then we could consider removing it from Octave as well, but we might want to warn that the feature is deprecated for a few releases before removing as we try to do for other interface changes.  FWIW, I was never a big fan of passing extra parameters this way and it wouldn't bother me personally to just remove the feature but I also know that users don't like backwardly incompatible changes.

Do you have a patch (even if preliminary) for this change in Octave?

John W. Eaton <jwe>
Group administrator
Wed 05 Oct 2022 05:29:19 PM UTC, comment #1: 

The solution to this bug involves sending additional arguments to file ode_event_handler.m so that it can call the interpolation function runge_kutta_interpolate.m. While I was making these changes, I fixed a couple other problems related to function parameters being passed around:

  • Extra function arguments were being passed around the ode solvers and sent to the derivative function (fcn), as well as OutputFcn and the Events function, if either of them existed. The extra arguments feature is undocumented and does not appear to exist in Matlab, and I believe the user is able to wrap their functions in anonymous functions to achieve the same result if desired. I am removing the sending of extra arguments to OutputFcn and Events, as it seems like poor practice to always send extra fcn arguments to other functions as well (especially when those functions are not documented to accept the extra arguments either). I am leaving the extra arguments feature for fcn at present, since it doesn't really seem to break any other functionality that way, although I have added a FIXME mark indicating it should possibly be removed.


  • ode_event_handler.m specifies that input y can be either a column vector, or a cell array with the y values in the first cell, and derivatives in the second. I could not find any function that calls ode_event_handler and utilizes the cell input. Interpolation would require handling this input if given, but since no functions actually use the cell input, I took it out of ode_event_handler and didn't try to implement it for the interpolation, either.
Ken Marek <m5318>
Wed 05 Oct 2022 05:07:44 PM UTC, original submission:  

The "Events" option in ode45 etc. currently uses linear interpolation to find Event locations (t and y), even though there is a Runge-Kutta interpolation function (of the appropriate polynomial order for the chosen solver) ready to use. This results in inaccurate location of the Events.

I have attached a zip file to demonstrate this behavior and compare to Matlab. test_file.m runs ode45 for two polynomial functions with known analytic solutions, order 2 and 7. Since ode45 uses a 4th order interpolation polynomial, we expect very good results for the order 2 problem, and some error in the order 7 problem.


Matlab results:
>> test_file
Elapsed time is 0.078885 seconds. % cleared memory
Elapsed time is 0.014720 seconds. % subsequent run
Elapsed time is 0.002338 seconds. % later runs
y2 function errors:
t error of event: 2.2204e-16
y error of "true" root: 4.4409e-16
y error of Event root, as reported by event: 1.1102e-16
y error of Event root, from y2(t_event): 6.6613e-16
y7 function errors:
t error of event: 2.6385e-07  6.3041e-08
y error of "true" root: 3.1086e-15  4.2633e-14
y error of Event root, as reported by event: 1.1102e-16  9.2015e-13
y error of Event root, from y7(t_event): 1.4023e-06  2.2269e-06



Octave results:
octave:35> test_file
Elapsed time is 0.0752091 seconds. % cleared memory
Elapsed time is 0.05669 seconds. % subsequent run
y2 function errors:
t error of event: 0.0067247
y error of "true" root: 4.4409e-16
y error of Event root, as reported by event: 1.1102e-16
y error of Event root, from y2(t_event): 0.035926
y7 function errors:
t error of event: 0.0063352    0.013498
y error of "true" root: 3.1086e-15  4.2633e-14
y error of Event root, as reported by event: 1.1102e-16  3.5527e-15
y error of Event root, from y7(t_event): 0.033534     0.46527


The Octave results do not have an additional time improvement for the third or later runs. (I am running the Matlab GUI and Octave on command line, both on Windows 10.)

There are also graphs at the same zoom level for the single root of the "y2" function, and for the two roots of the "y7" function, to visually compare Matlab and Octave results. We see that the true root is basically on the solution curve produced by ode45 for both programs, but Octave's Event handler adds significantly more error than Matlab's. (You may notice in the Octave plots that the Event result is on the other side of a solution point from the true solution, which might seem problematic. However, the point shown is an interpolation "Refine" point, not a base solution point, so the Event result is actually bounded by the same two base solution points as the true result.) You can also see this in the text output; the most relevant error measurement is the "y error of Event root, from yn(t_event)". We can see that for the "y2" function, Octave returns y results that are 14 orders of magnitude worse than Matlab; for "y7", they are "only" 5 orders worse. Results for t are similar. The problem resides mostly in file ode_event_handler.m, so it affects ode23 and ode23s as well (I have verified this but am not posting additional details; you can change the solver in the test file to verify).

I have a solution ready to post, but I am waiting on the solution to bug #63063 to be pushed. This touches a lot of the same files, and would possibly break without that patch being applied first. My patch will also clean up some other outdated or unused parts of the code, which I will describe in another comment.

Ken Marek <m5318>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #53811:  event_location_bug63162_20221006.patch added by m5318 (13KiB - application/octet-stream)
file #53810:  rkinterp_bug63163_20221006.patch added by m5318 (3KiB - application/octet-stream)
file #53797:  patch_bug63162_20221005.patch added by m5318 (17KiB - application/octet-stream)
file #53796:  event_interpolation.zip added by m5318 (99KiB - application/x-zip-compressed)

 

Digest:
   bug dependencies.

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by cdf (Posted a comment)
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by m5318 (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 11 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-11-21 nrjank StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2022-11-11 nrjank StatusNeed Info Ready For Test
    2022-11-10 nrjank StatusIn Progress Need Info
    2022-10-06 m5318 Attached File- Added event_location_bug63162_20221006.patch, #53811
    2022-10-06 m5318 Attached File- Added rkinterp_bug63163_20221006.patch, #53810
    2022-10-05 nrjank StatusNone In Progress
        Release7.2.0 dev
        Dependencies- Depends on bugs #63063
    2022-10-05 m5318 Attached File- Added patch_bug63162_20221005.patch, #53797
    2022-10-05 m5318 Attached File- Added event_interpolation.zip, #53796

    Back to the top

    Powered by Savane 3.12.
    Corresponding source code