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/api2.2/org/apache/commons/math/analysis/solvers/packagesummary.html.
>> In particular, the UnivariateRealSolver interface seems to be the
>> interface for all rootfinding 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
>> rootfinding algorithm. It seems that the BrentSolver rootfinding
>> algorithm is hardcoded, and no other rootfinding algorithm can be used
>> to find state events...
>
> Yes, it is hardcoded.
>> I was wondering if there is any possibility to use another rootfinding
>> 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 GraggBulirschStoer, 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(1e10, 100.0, 1e10, 1e10)
>>>> integrator.addEventHandler(event, 100.0, 1e10, 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, email: userunsubscribe@commons.apache.org
>>>> For additional commands, email: userhelp@commons.apache.org
>>> 
>>> To unsubscribe, email: userunsubscribe@commons.apache.org
>>> For additional commands, email: userhelp@commons.apache.org
>>>
>>
>> 
>> To unsubscribe, email: userunsubscribe@commons.apache.org
>> For additional commands, email: userhelp@commons.apache.org
>>
>
>
> 
> To unsubscribe, email: userunsubscribe@commons.apache.org
> For additional commands, email: userhelp@commons.apache.org
>

To unsubscribe, email: userunsubscribe@commons.apache.org
For additional commands, email: userhelp@commons.apache.org
