[ https://issues.apache.org/jira/browse/MATH705?page=com.atlassian.jira.plugin.system.issuetabpanels:commenttabpanel&focusedCommentId=13147240#comment13147240
]
Luc Maisonobe edited comment on MATH705 at 11/9/11 7:58 PM:

The reason for this strange behavior is that g function evaluations are based on the integratorspecific
interpolator.
Each integration method has its specific algorithm and preserve a rich internal data set.
From this data set, it is possible to build an interpolator which is specific to the integrator
and in fact share part of the data set (they reference the same arrays). So integrator and
interpolator are tightly bound together.
For embedded RungeKutta methods like DormandPrince 8(5,3), this data set corresponds to
one state vector value and several state vector derivatives sampled throughout the step. When
the step is accepted after error estimation, the state value is set to the value at end of
step and the interpolator is called. So the equations of the interpolator are written in such
a way interpolation is backward: we start from the end state and roll back to beginning of
step. This explains why when we roll all the way back to step start, we may find a state that
is not exactly the one we started from, due to both the integration order and interpolation
order.
For GraggBulirschStoer, the data set corresponds to one state vector value and derivatives
at several orders, all taken at step middle point. When the step is accepted after error estimation,
the interpolator is called before the state value is set to the value at end of step and.
So the equations of the interpolator are written in such a way interpolation is forward: we
start from the start state and go on towards end of step. This explains why when we go all
the way to step end, we may find a state that is not exactly the one that will be used for
next step, due to both the integration order and interpolation order.
So one integrator type is more consistent at step start and has more error at step end, while
the other integrator has a reversed behavior.
In any case, the interpolation that is used (and in fact the integration data set it is based
upon) are not error free. The error is related to step size.
We could perhaps rewrite some interpolators by preserving both start state s(t[k]) and end
state s(t[k+1]) and switching between two hal model as follows:
i(t) = s(t[k]) + forwardModel(t[k], t) if t <= (t[k] + t[k+1])/2
and
i(t) = s(t[k+1]) + backwardModel(t, t[k+1]) if t > (t[k] + t[k+1])/2
This would make interpolator more consistent with integrator at both step start and step end
and perhaps reduce this problem. This would however not be perfect, as it will introduce a
small error at junction point. I'm not sure if it would be easy or not, we would have to review
all interpolators and all integrators for that. All models are polynomial ones.
Note that the problem should not appear for Adams methods (when they will be considered validated
...), because in this case, it is the interpolator that is built first and the integrator
is in fact an application of the interpolator at step end! So interpolator and integrator
are by definition always perfectly consistent with each other.
What do you think ?
Should we let this problem alone and consider we are in the grey zone of expected numerical
inaccuracy due to integration/interpolation orders or should we attempt the two halfmodels
trick ?
was (Author: luc):
The reason for this strange behavior is that g function evaluations are based on the integratorspecific
interpolator.
Each integration method has its specific algorithm and preserve a rich internal data set.
From this data set, it is possible to build an interpolator which is specific to the integrator
and in fact share part of the data set (they reference the same arrays). So integrator and
interpolator are tightly bound together.
For embedded RungeKutta methods like DormandPrint 8(5,3), this data set corresponds to one
state vector value and several state vector derivatives sampled throughout the step. When
the step is accepted after error estimation, the state value is set to the value at end of
step and the interpolator is called. So the equations of the interpolator are written in such
a way interpolation is backward: we start from the end state and roll back to beginning of
step. This explains why when we roll all the way back to step start, we may find a state that
is not exactly the one we started from, due to both the integration order and interpolation
order.
For GraggBulirschStoer, the data set corresponds to one state vector value and derivatives
at several orders, all taken at step middle point. When the step is accepted after error estimation,
the interpolator is called before the state value is set to the value at end of step and.
So the equations of the interpolator are written in such a way interpolation is forward: we
start from the start state and go on towards end of step. This explains why when we go all
the way to step end, we may find a state that is not exactly the one that will be used for
next step, due to both the integration order and interpolation order.
So one integrator type is more consistent at step start and has more error at step end, while
the other integrator has a reversed behavior.
In any case, the interpolation that is used (and in fact the integration data set it is based
upon) are not error free. The error is related to step size.
We could perhaps rewrite some interpolators by preserving both start state s(t[k]) and end
state s(t[k+1]) and switching between two hal model as follows:
i(t) = s(t[k]) + forwardModel(t[k], t) if t <= (t[k] + t[k+1])/2
and
i(t) = s(t[k+1]) + backwardModel(t, t[k+1]) if t > (t[k] + t[k+1])/2
This would make interpolator more consistent with integrator at both step start and step end
and perhaps reduce this problem. This would however not be perfect, as it will introduce a
small error at junction point. I'm not sure if it would be easy or not, we would have to review
all interpolators and all integrators for that. All models are polynomial ones.
Note that the problem should not appear for Adams methods (when they will be considered validated
...), because in this case, it is the interpolator that is built first and the integrator
is in fact an application of the interpolator at step end! So interpolator and integrator
are by definition always perfectly consistent with each other.
What do you think ?
Should we let this problem alone and consider we are in the grey zone of expected numerical
inaccuracy due to integration/interpolation orders or should we attempt the two halfmodels
trick ?
> DormandPrince853 integrator leads to revisiting of state events
> 
>
> Key: MATH705
> URL: https://issues.apache.org/jira/browse/MATH705
> Project: Commons Math
> Issue Type: Bug
> Affects Versions: 3.0
> Environment: Commons Math trunk, Java 6, Linux
> Reporter: Dennis Hendriks
> Attachments: ReappearingEventTest.java, ReappearingEventTest.out
>
>
> See the attached ReappearingEventTest.java. It has two unit tests, which use either the
DormandPrince853 or the GraggBulirschStoer integrator, on the same ODE problem. It is a problem
starting at time 6.0, with 7 variables, and 1 state event. The state event was previously
detected at time 6.0, which is why I start there now. I provide and end time of 10.0. Since
I start at the state event, I expect to integrate all the way to the end (10.0). For the GraggBulirschStoer
this is what happens (see attached ReappearingEventTest.out). For the DormandPrince853Integerator,
it detects a state event and stops integration at 6.000000000000002.
> I think the problem becomes clear by looking at the output in ReappearingEventTest.out,
in particular these lines:
> {noformat}
> computeDerivatives: t=6.0 y=[2.0 , 2.0
, 2.0 , 4.0 , 2.0 , 7.0
, 15.0 ]
> (...)
> g : t=6.0 y=[1.9999999999999996 , 1.9999999999999996
, 1.9999999999999996 , 4.0 , 1.9999999999999996 , 7.0 ,
14.999999999999998 ]
> (...)
> final result : t=6.000000000000002 y=[2.0000000000000013 , 2.0000000000000013
, 2.0000000000000013 , 4.000000000000002 , 2.0000000000000013 , 7.000000000000002 ,
15.0 ]
> {noformat}
> The initial value of the last variable in y, the one that the state event refers to,
is 15.0. However, the first time it is given to the g function, the value is 14.999999999999998.
This value is less than 15, and more importantly, it is a value from the past (as all functions
are increasing), *before* the state event. This makes that the state event reappears immediately,
and integration stops at 6.000000000000002 because of the detected state event.
> I find it puzzling that for the DormandPrince853Integerator the y array that is given
to the first evaluation of the g function, has different values than the y array that is the
input to the problem. For GraggBulirschStoer is can be seen that the y arrays have identical
values.

This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira
