[Octave-bug-tracker] [bug #63162] Inaccurate "Events" location in ode45, ode23(s)

[Ken Marek]

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

                 Summary: Inaccurate "Events" location in ode45, ode23(s)
                 Project: GNU Octave
               Submitter: 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: None
             Assigned to: None
         Originator Name: 
        Originator Email: 
             Open/Closed: Open
                 Release: 7.2.0
         Discussion Lock: Any
        Operating System: Any


    _______________________________________________________

Follow-up Comments:


-------------------------------------------------------
Date: Wed 05 Oct 2022 05:07:44 PM UTC By: Ken Marek <m5318>
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.






    _______________________________________________________
File Attachments:


-------------------------------------------------------
Date: Wed 05 Oct 2022 05:07:44 PM UTC  Name: event_interpolation.zip  Size:
99KiB   By: m5318

<http://savannah.gnu.org/bugs/download.php?file_id=53796>

    _______________________________________________________

Reply to this item at:

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

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