bugGNU Octave - Bugs: bug #49919, ode45 event problem, event active...

 
 

bug #49919: ode45 event problem, event active at initial time

Submitter:  Bill Greene <billgreene>
Submitted:  Fri 23 Dec 2016 02:50:08 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Matlab Compatibility
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

Sun 22 Jan 2017 12:43:41 AM UTC, comment #9: 

I checked in a fix here (http://hg.savannah.gnu.org/hgweb/octave/rev/17dc6c7ef427).  This was not easy, although the patch looks easy in the end.  The regression tests are passing, as well as the two tests in this bug report, but feel free to do more testing.

Rik <rik5>
Group administrator
Wed 04 Jan 2017 08:11:06 PM UTC, comment #8: 

r2017a prerelease:

>> event_first_step2
te =
    0.0000
    2.3509
ye =
    1.0000    0.0000   -0.0000    0.3000
    0.9811   -0.0001   -0.0020    0.3060
ie =
     1
     1
tfinal =
    2.3509
yfinal =
    0.9811   -0.0001   -0.0020    0.3060
>>


+ attached figure.


Philip Nienhuis <philipnienhuis>
Group Member
Wed 04 Jan 2017 06:01:23 PM UTC, comment #7: 

Does anyone have access to Matlab in order to run the test file in comment #6?  I need that to make progress.

Rik <rik5>
Group administrator
Mon 26 Dec 2016 07:02:26 AM UTC, comment #6: 

This turns out to be quite an involved problem.  In order to determine exactly how Matlab is handling the event function, could someone with access to Matlab run the attached event_first_step2.m file and report the results?


(file #39313)

Rik <rik5>
Group administrator
Fri 23 Dec 2016 07:08:08 PM UTC, comment #5: 

I just realized that I mis-labeled Release when I filed this bug--
sorry. I am running my tests with development code after your change for
 bug #49846.

I get

tfinal = 0
yfinal =

   1.00000   0.00000   0.00000   0.30000

which is also wrong. This result shows the incorrect termination
at t=0 and yfinal = the initial conditions.

Bill Greene <billgreene>
Fri 23 Dec 2016 06:03:42 PM UTC, comment #4: 

Possibly I have mis-diagnosed the root cause of the problem.
But, yes, I think bug #49846 is either the same problem or
closely related.

A matlab doc that discusses this "first step" issue is
https://www.mathworks.com/help/matlab/math/ode-event-location.html
under the section "Limitations"

This example is from Moler's book Numerical Computing With Matlab
If you look on page 26 of this doc
https://www.mathworks.com/content/dam/mathworks/mathworks-dot-com/moler/odes.pdf
you will see what (apparently) matlab returns for this example.
As you can see from the figure, matlab terminates when the solution
gets back to y=0 the second time around.

Bill Greene <billgreene>
Fri 23 Dec 2016 05:53:26 PM UTC, comment #3: 

Okay, the issue seems to be the second one.  See scritps/ode/private/integrate_adaptive.m.  The function ode_event_handler is called outside of the main integration to initialize everything.  And then the main loop is entered and the function is called again at time0.

Rik <rik5>
Group administrator
Fri 23 Dec 2016 05:50:54 PM UTC, comment #2: 

The function in question is in scripts/ode/private/ode_event_handler.m.  It is supposed to be called at the start of the integration with the 'init' argument so it can initialize the internal variables.

Either we are not calling it with the 'init' option, or we call it once before the integration starts with the 'init' option and then call it again, without 'init', at time 0 in the integration.

Rik <rik5>
Group administrator
Fri 23 Dec 2016 05:46:01 PM UTC, comment #1: 

@Bill: Can you send a pointer to the location in the Matlab docs where this is explained?  I don't think it was clear what we were aiming for.

Also, I think bug #49846 (wrong initial time for Event function) bears on this.  If I run your test script with Octave-4.2.0 I get


event_first_step
ie =  1
tfinal =  6.2832
yfinal =

   1.0000e+00   2.7756e-17  -1.1102e-16   3.0000e-01


If I run the stable branch code after


changeset:   22894:544dac72d414
branch:      stable
parent:      22891:fd056e736f73
user:        Rik <rik@octave.org>
date:        Thu Dec 15 15:23:43 2016 -0800
summary:     Initialize ODE Event function with start time (bug #49846).


I get


octave:1> event_first_step
ie =  1
tfinal = 0
yfinal =

   1.0000e+00   2.7756e-17  -1.1102e-16   3.0000e-01


The variable tfinal looks correct because the last integration time was 0.  I'm tempted to close #49846.


Rik <rik5>
Group administrator
Fri 23 Dec 2016 02:50:08 PM UTC, original submission:  

This is a bug I've been working on but don't yet have a solution.
So I'm submitting this hoping the solution is obvious to someone
with more experience with this code.

In this example, the event condition is active at the initial time and the code terminates in the first step. But according to matlab docs, the first step is supposed to be handled specially so that termination doesn't occur until the next time the event becomes active.

I've attached a simple example of this.

Bill Greene <billgreene>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #39313:  event_first_step2.m added by rik5 (706B - d2l/unknowntype)
file #39300:  event_first_step.m added by billgreene (626B - text/plain)

 

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 billgreene (Submitted the item)
  • -email is unavailable- added by billgreene
  •  

    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
    2017-01-22 rik5 StatusNeed Info Fixed
        Open/ClosedOpen Closed
    2017-01-04 philipnienhuis Attached File- Added event_first_step2.png, #39365
    2016-12-30 rik5 StatusConfirmed Need Info
    2016-12-26 rik5 Attached File- Added event_first_step2.m, #39313
    2016-12-23 rik5 Item GroupIncorrect Result Matlab Compatibility
        StatusNeed Info Confirmed
        Release4.2.0 dev
    2016-12-23 rik5 StatusNone Need Info
    2016-12-23 billgreene Attached File- Added event_first_step.m, #39300
        Carbon-Copy- Added -email is unavailable-

    Back to the top

    Powered by Savane 3.13-caa5.
    Corresponding source code