----- "Dennis Hendriks" <d.hendr...@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-unsubscr...@commons.apache.org
For additional commands, e-mail: user-h...@commons.apache.org
---------------------------------------------------------------------
To unsubscribe, e-mail: user-unsubscr...@commons.apache.org
For additional commands, e-mail: user-h...@commons.apache.org