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 08:47:21 GMT
Hi Luc,

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?

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


Mime
View raw message