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 Fri, 10 Jun 2011 10:50:46 GMT
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


Mime
View raw message