[jira] [Issue Comment Edited] (MATH-705) DormandPrince853 integrator leads to revisiting of state events
[ https://issues.apache.org/jira/browse/MATH-705?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanelfocusedCommentId=13147712#comment-13147712 ] Dennis Hendriks edited comment on MATH-705 at 11/10/11 2:31 PM: bq. 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 half-models trick? I consider it important that events are properly detected, and are detected exactly once (if they occur exactly once). Therefore, in general, I think the half-models trick would be a good idea, as it is more important (to me) to have higher accuracy at the end points than in the middle. For me, the DormandPrince853Integrator is now practically useless. was (Author: dhendriks): 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 half-models trick? I consider it important that events are properly detected, and are detected exactly once (if they occur exactly once). Therefore, in general, I think the half-models trick would be a good idea, as it is more important (to me) to have higher accuracy at the end points than in the middle. For me, the DormandPrince853Integrator is now practically useless. DormandPrince853 integrator leads to revisiting of state events --- Key: MATH-705 URL: https://issues.apache.org/jira/browse/MATH-705 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.002. 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.9996 , 1.9996 , 1.9996 , 4.0 , 1.9996 , 7.0 , 14.998 ] (...) final result : t=6.002y=[2.0013 , 2.0013 , 2.0013 , 4.002 , 2.0013 , 7.002 , 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.998. 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 re-appears immediately, and integration stops at 6.002 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
[jira] [Issue Comment Edited] (MATH-705) DormandPrince853 integrator leads to revisiting of state events
[ https://issues.apache.org/jira/browse/MATH-705?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanelfocusedCommentId=13147240#comment-13147240 ] Luc Maisonobe edited comment on MATH-705 at 11/9/11 7:58 PM: - The reason for this strange behavior is that g function evaluations are based on the integrator-specific 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 Runge-Kutta methods like Dormand-Prince 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 Gragg-Bulirsch-Stoer, 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 half-models trick ? was (Author: luc): The reason for this strange behavior is that g function evaluations are based on the integrator-specific 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 Runge-Kutta methods like Dormand-Print 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 Gragg-Bulirsch-Stoer, 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