commons-user mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Dennis Hendriks <D.Hendr...@tue.nl>
Subject Re: [Math] How to force integration past the state event?
Date Wed, 15 Jun 2011 07:15:53 GMT
Just to let everyone know: the issue 
(https://issues.apache.org/jira/browse/MATH-586) is now closed/fixed, and a 
patch was committed to trunk as of r1135726.

Best regards,
Dennis


Dennis Hendriks wrote:
> Hi Luc,
> 
> I just created an issue and attached a patch to it. See 
> https://issues.apache.org/jira/browse/MATH-586
> 
> Is there anything else I should do?
> 
> Best regards,
> Dennis
> 
> 
> Luc Maisonobe wrote:
>> Le 10/06/2011 10:47, Dennis Hendriks a écrit :
>>> Hi Luc,
>> Hi Dennis,
>>
>>> Thanks again. I'm willing to try an make a patch for this. I currently
>>> used release 2.2. I assume I would have to switch to the current
>>> development version that is to become 3.0?
>> Yes. This change would imply modifying an existing top level public 
>> interface to add the set method, so it can be done only when a major 
>> release is published. So it is right the perfect time to do it!
>>
>> best regards,
>> Luc
>>
>>> Best regards,
>>> Dennis
>>>
>>>
>>> Luc Maisonobe wrote:
>>>> Le 10/06/2011 10:20, Dennis Hendriks a écrit :
>>>>> Hi Luc,
>>>> Hi Dennis,
>>>>
>>>>> Thanks for your quick reply. It would be a possible solution.
>>>>>
>>>>> I decided to familiarize myself with the internals of the code that is
>>>>> used, in order to better understand what code is responsible for the
>>>>> root finding. I found the org.apache.commons.math.analysis.solvers
>>>>> package, at
>>>>> http://commons.apache.org/math/api-2.2/org/apache/commons/math/analysis/solvers/package-summary.html.
>>>>>
>>>>> In particular, the UnivariateRealSolver interface seems to be the
>>>>> interface for all root-finding algorithms.
>>>> Yes.
>>>>
>>>>> I also found the EventState class
>>>>> (http://svn.apache.org/viewvc/commons/proper/math/tags/MATH_2_2/src/main/java/org/apache/commons/math/ode/events/EventState.java?view=markup)
>>>>>
>>>>> which at line 242 creates an instance of the BrentSolver to use as
>>>>> root-finding algorithm. It seems that the BrentSolver root-finding
>>>>> algorithm is hard-coded, and no other root-finding algorithm can be used
>>>>> to find state events...
>>>> Yes, it is hard-coded.
>>>>> I was wondering if there is any possibility to use another root-finding
>>>>> algorithm for state event detection during ODE solver integration. That
>>>> Allowing to customize this would be an improvement. This should be
>>>> done at integrator level probably using a setEventRootSolver method
>>>> for example and passed down to EventState. The EventState class by
>>>> itself is not managed by users. Also this method should be implemented
>>>> for all integrators types, not only Gragg-Bulirsch-Stoer, but it is
>>>> quite simple.
>>>>
>>>>> would allow me to create my own UnivariateRealSolver implementation that
>>>>> does guarantee that the point returned never undershoots the state
>>>>> event. I believe that would be the preferred solution to my initial
>>>>> problem.
>>>> Yes, it is a better solution. If you think you could contribute the
>>>> library change, do not hesitate to open an issue in our tracker and
>>>> attache a patch to it.
>>>>
>>>> Thanks,
>>>> Luc
>>>>
>>>>> Best regards,
>>>>> Dennis
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> luc.maisonobe@free.fr wrote:
>>>>>> ----- "Dennis Hendriks" <D.Hendriks@tue.nl> a écrit :
>>>>>>
>>>>>>> Hi all,
>>>>>> Hi Dennis,
>>>>>>
>>>>>>> I have the following (simplified) ODE problem:
>>>>>>>
>>>>>>> V(0) = 10.000000645817822
>>>>>>> V' = -sqrt(V)
>>>>>>>
>>>>>>> and the following state event:
>>>>>>>
>>>>>>> V <= 2.0
>>>>>>>
>>>>>>> Using:
>>>>>>>
>>>>>>> t = 6.6845103160078425
>>>>>>> GraggBulirschStoerIntegrator(1e-10, 100.0, 1e-10, 1e-10)
>>>>>>> integrator.addEventHandler(event, 100.0, 1e-10, 999)
>>>>>>>
>>>>>>> I obtain the last time point as:
>>>>>>>
>>>>>>> t = 10.180638128882343
>>>>>>> V = 2.000000830098894
>>>>>>> V' = -1.414213855857346
>>>>>>>
>>>>>>> For this last time point, V <= 2.0 does not hold. In my previous
>>>>>> What you experience here is that convergence is reached on the side
of
>>>>>> the root you don't want.
>>>>>> From a pure event root finder point of view, there is no guarantee
>>>>>> about the
>>>>>> side at which the solver will stop. You can reduce the convergence
>>>>>> threshold, of course,
>>>>>> but it will only make the undershoot smaller, it will not always
>>>>>> prevent it.
>>>>>>
>>>>>>> experiences with DDASRT (DAE solver written in Fortran, see
>>>>>>> http://www.netlib.org/ode/ddasrt.f), it was actually guaranteed
that
>>>>>>> integration would go PAST the state event.
>>>>>>>
>>>>>>> So, my question is: Is there any way to force integration PAST
the
>>>>>>> state event, to make sure that afterwards the V <= 2.0 condition
>>>>>>> actually
>>>>>>> holds? In other words, is it possible to make sure it stops PAST
the
>>>>>>> state
>>>>>>> event, instead of sometimes PAST it and sometimes BEFORE it?
>>>>>> I would suggest to use a StepHandler in addition to the EventHandler,
>>>>>> and to
>>>>>> retrieve the final result from the step handler instead of the ODE
>>>>>> itself.
>>>>>> This means, instead of:
>>>>>>
>>>>>> integrator.addEventHandler(...);
>>>>>> double tEnd = integrator.integrate(...);
>>>>>>
>>>>>> you would use:
>>>>>>
>>>>>> integrator.addEventHandler(...);
>>>>>> integrator.addStepHandler(mySpecialStepHandler);
>>>>>> integrator.integrate(...);
>>>>>> double tEnd = mySpecialStepHandler.getTend();
>>>>>> double[] yEnd = mySpecialStepHandler.getY();
>>>>>>
>>>>>> with myStepHandler being a local class of your own with an
>>>>>> implementation like:
>>>>>>
>>>>>> public class MyStepHandler implements StepHandler {
>>>>>>
>>>>>> double tEnd;
>>>>>> double[] y;
>>>>>>
>>>>>> public boolean requiresDenseOutput() {
>>>>>> return true;
>>>>>> }
>>>>>>
>>>>>> public void reset() {
>>>>>> }
>>>>>>
>>>>>> public void handleStep(StepInterpolator interpolator, boolean isLast)
{
>>>>>> if (isLast) {
>>>>>> tEnd = interpolator.getCurrentTime();
>>>>>> while (true) {
>>>>>> interpolator.setInterpolatedTime(tEnd);
>>>>>> yEnd = interpolator.getInterpolatedState();
>>>>>> if (yEnd[0] <= 2) {
>>>>>> return;
>>>>>> }
>>>>>>
>>>>>> // we are on the wrong side, slightly shift end time
>>>>>> double[] yEndDot = interpolator.getInterpolatedDerivatives();
>>>>>> tEnd += (yEnd[0] - 2) / yEndDot[0];
>>>>>>
>>>>>> }
>>>>>> }
>>>>>> }
>>>>>>
>>>>>> public double getTend() {
>>>>>> return tEnd;
>>>>>> }
>>>>>>
>>>>>> public double[] getYend() {
>>>>>> return yEnd;
>>>>>> }
>>>>>>
>>>>>> }
>>>>>>
>>>>>> The trick here is that the interpolator can be used slightly out
of
>>>>>> the current step,
>>>>>> so slightly shifting the end date after it has been computed of the
>>>>>> order of magnitude
>>>>>> of the convergence threshold in the event handler is safe.
>>>>>>
>>>>>> It is guaranteed that the event handler is called after the step
>>>>>> handler, in fact it is
>>>>>> exactly because we need to make sure the isLast boolean is properly
>>>>>> set to true when an
>>>>>> events ask the integrator to stop.
>>>>>>
>>>>>>> Any help would be greatly appreciated.
>>>>>> Hope this helps,
>>>>>>
>>>>>> Luc
>>>>>>
>>>>>>> Thanks,
>>>>>>> Dennis
>>>>>>>
>>>>>>> ---------------------------------------------------------------------
>>>>>>> To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
>>>>>>> For additional commands, e-mail: user-help@commons.apache.org
>>>>>> ---------------------------------------------------------------------
>>>>>> To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
>>>>>> For additional commands, e-mail: user-help@commons.apache.org
>>>>>>
>>>>> ---------------------------------------------------------------------
>>>>> To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
>>>>> For additional commands, e-mail: user-help@commons.apache.org
>>>>>
>>>> ---------------------------------------------------------------------
>>>> To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
>>>> For additional commands, e-mail: user-help@commons.apache.org
>>>>
>>> ---------------------------------------------------------------------
>>> To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
>>> For additional commands, e-mail: user-help@commons.apache.org
>>>
>>
>> ---------------------------------------------------------------------
>> To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
>> For additional commands, e-mail: user-help@commons.apache.org
>>
> 
> 
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
> For additional commands, e-mail: user-help@commons.apache.org
> 


---------------------------------------------------------------------
To unsubscribe, e-mail: user-unsubscribe@commons.apache.org
For additional commands, e-mail: user-help@commons.apache.org


Mime
View raw message