[Octave-bug-tracker] [bug #63063] ode{23, 45} options "Events" and "OutputFcn" do not perform correctly.

[Ken Marek]

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

                 Summary: ode{23,45} options "Events" and "OutputFcn" do not
perform correctly.
                 Project: GNU Octave
               Submitter: m5318
               Submitted: Thu 15 Sep 2022 09:08:18 PM UTC
                Category: Octave Function
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Incorrect Result
                  Status: None
             Assigned to: None
         Originator Name: 
        Originator Email: 
             Open/Closed: Open
                 Release: 7.2.0
         Discussion Lock: Any
        Operating System: Microsoft Windows


    _______________________________________________________

Follow-up Comments:


-------------------------------------------------------
Date: Thu 15 Sep 2022 09:08:18 PM UTC By: Ken Marek <m5318>
Overview: Each of these options displays incorrect behavior on its own.
Additionally, there is an interaction between the two bugs, so I am treating
this as a single bug. (Interaction is: termination due to an Event prevents
OutputFcn from being called on the last iteration.) I am currently working on
a fix to bug #49408 relating to ode{23,45} "Refine", and that fix will require
fixing these issues as well. This bug report is to document these additional
problems.

To display these bugs, we need two functions defined. I have made the
following simple examples:

function [val, isTerm, dir] = testEvent ( t , y )
  isTerm = [ 0; 1 ];
  dir = [ 1; 1 ];
  val = [ y(1) - 1.1; y(1) - 1.5 ];
end


and


function stopSolve = testFuncOut ( t, y, flag )
  persistent stepNum  
  stopVal = 1.2;
  if strcmp ( flag, 'init' )
    stepNum = 0;
  elseif strcmp ( flag, 'done' )
    disp ( 'Output function received DONE signal' )
    return
  else
    stepNum = stepNum + 1;
  end
  disp ( ['Step: ' num2str(stepNum) ', t = ' num2str(t(1)) ', y = ' num2str(y)
] )
  if y(1) > stopVal
    stopSolve = true;
    disp ( ['y value is ' num2str(y(1)) ', greater than ' num2str(stopVal) '.
Sending STOP signal.'] )
  else
    stopSolve = false;
  end
end


These examples are in Matlab compatible script so we can compare results.

Function testEvent creates a non-terminating event at y = 1.1, and a
terminating event at y = 1.5 (for scalar y). Function testFuncOut prints t and
y to the screen each time the function is called, except when the "done" flag
is sent, in which case it indicates it received the "done" signal.
Additionally, testFuncOut returns a termination signal if y > 1.2.

Now I have created 4 test cases:

%test 1 (baseline)
odeOpt = odeset ( 'Events', [], 'InitialStep', 1, 'OutputFcn', [], ...
  'Refine', 1, 'MaxStep', 1 );

%test 2 (event only)
odeOpt = odeset ( 'Events', @testEvent, 'InitialStep', 1, 'OutputFcn', [],
...
  'Refine', 1, 'MaxStep', 1 );

%test 3 (output function only)
odeOpt = odeset ( 'Events', [], 'InitialStep', 1, 'OutputFcn', @testFuncOut,
...
  'Refine', 1, 'MaxStep', 1 );

%test 4 (event and output function)
odeOpt = odeset ( 'Events', @testEvent, 'InitialStep', 1, 'OutputFcn',
@testFuncOut, ...
  'Refine', 1, 'MaxStep', 1 );


With each set of odeOpt options, I ran ode45 using:

[t, y] = ode45 ( @(t,y) 1, [0,1,2,2.5], [0], odeOpt )

I also ran with struct output:

structOut = ode45 ( @(t,y) 1, [0,1,2,2.5], [0], odeOpt )

when needed to pull out the returned Event variables (one could alternately
assign additional return variables to get them). The ODE I have selected has
the solution y(t) = t, so ode45 will return only exact solution results.

Results are:
----------
Test 1

Matlab:
t = [0; 1; 2; 2.5]
y = [0; 1; 2; 2.5]
Octave: same

There are no errors for the test scenario when Events and OutputFcn options
are not used.

----------
Test 2
Matlab:
t = [0; 1; 1.5]
y = [0; 1; 1.5]
event variables:
  xe = [1.1, 1.5]
  ye = [1.1, 1.5]
  ie = [1, 2]
Octave: produces warning. Also:
t = [0; 1; 1.1]
y = [0; 1; 1.1]
event variables:
  xe = [1.1]
  ye = [1.1]
  ie = [1]

Here we see an error in the Events option alone. Octave produces a warning
where Matlab does not, and it also incorrectly treats the non-terminal Event
as a terminal one. I do not plan to remove the warning in my bug fixes, unless
asked to do so.

----------
Test 3
Matlab:
t = [0; 1; 2]
y = [0; 1; 2]
testFuncOut shows steps 0, 1, and 2
Octave: produces warning, otherwise same

There is no error here when OutputFcn is used alone (but we will see an
OutputFcn bug in a later test case).

----------
Test 4
Matlab:
t = [0; 1; 1.5]
y = [0; 1; 1.5]
testFuncOut shows steps 0, 1, and 2
Octave: produces warning. Also:
t = [0; 1; 1.1]
y = [0; 1; 1.1]
testFuncOut shows only steps 0 and 1

Here, termination due to an Event prevents OutputFcn from being called on the
last step in Octave.


----------

The above function calls to ode45 used a fixed set of output times. We now try
again using a 2-element trange:

[t, y] = ode45 ( @(t,y) 1, [0,2.5], [0], odeOpt )

Because of the initial step and maximum step size I have chosen, the results
here should have the same t values as before, so all results should
theoretically be identical to the above tests. Actual results are:

----------
Tests 1-4 for Matlab and Tests 1-2 for Octave are the same as above.

----------
Test 3 Octave:
testFuncOut shows steps numbered 0 to 4 (t = 0 and t = 1 are each repeated
once).

This is an error in OutputFcn with no Events, and is due to an incorrect
implementation of Refine (hence the connection to the other bug, as mentioned
in the summary).

----------
Test 4 Octave:
Results are qualitatively similar to Test 3. testFuncOut shows steps numbered
0 to 2 (t = 0 is repeated once). The Events/OutputFcn interaction bug is also
still present.

----------

A brief search did not reveal this to be an open bug, though I found several
bugs about Events initialization, which I think have all been fixed. I think
they are all in ode_event_handler.m, which I don't expect to modify in these
fixes, so hopefully I won't break anything in the process.

As far as I have identified, all the bugs are in the file
integrate_adaptive.m, which is called by both ode45 and ode23. Therefore, bugs
as shown here with ode45 will be also be expressed in ode23.








    _______________________________________________________

Reply to this item at:

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

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