Mailing-List: contact user-help@commons.apache.org; run by ezmlm
Precedence: bulk
Reply-To: "Commons Users List" <user@commons.apache.org>
Received-SPF: neutral (athena.apache.org: local policy)
Message-ID: <4DF1D5FF.7060402@free.fr>
Date: Fri, 10 Jun 2011 10:29:51 +0200
From: Luc Maisonobe <Luc.Maisonobe@free.fr>
User-Agent: Mozilla/5.0 (X11; U; Linux x86_64; en-US;
 rv:1.9.2.15) Gecko/20110419 Thunderbird/3.1.9
MIME-Version: 1.0
To: Commons Users List <user@commons.apache.org>
Subject: Re: [Math] How to force integration past the state event?
References: 
 <1046701180.1582351307636302599.JavaMail.root@spooler6-g27.priv.proxad.net>
 <4DF1D3D8.30004@tue.nl>
In-Reply-To: <4DF1D3D8.30004@tue.nl>
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit

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