mirror of
https://github.com/apache/commons-math.git
synced 2025-02-06 01:59:13 +00:00
added event handling to ODE integrators with jacobians
git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@919844 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
9034808411
commit
f5f8956448
@ -0,0 +1,203 @@
|
||||
/*
|
||||
* Licensed to the Apache Software Foundation (ASF) under one or more
|
||||
* contributor license agreements. See the NOTICE file distributed with
|
||||
* this work for additional information regarding copyright ownership.
|
||||
* The ASF licenses this file to You under the Apache License, Version 2.0
|
||||
* (the "License"); you may not use this file except in compliance with
|
||||
* the License. You may obtain a copy of the License at
|
||||
*
|
||||
* http://www.apache.org/licenses/LICENSE-2.0
|
||||
*
|
||||
* Unless required by applicable law or agreed to in writing, software
|
||||
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
* See the License for the specific language governing permissions and
|
||||
* limitations under the License.
|
||||
*/
|
||||
|
||||
package org.apache.commons.math.ode.jacobians;
|
||||
|
||||
import org.apache.commons.math.ode.events.EventException;
|
||||
|
||||
/** This interface represents a handler for discrete events triggered
|
||||
* during ODE integration.
|
||||
*
|
||||
* <p>Some events can be triggered at discrete times as an ODE problem
|
||||
* is solved. This occurs for example when the integration process
|
||||
* should be stopped as some state is reached (G-stop facility) when the
|
||||
* precise date is unknown a priori, or when the derivatives have
|
||||
* discontinuities, or simply when the user wants to monitor some
|
||||
* states boundaries crossings.
|
||||
* </p>
|
||||
*
|
||||
* <p>These events are defined as occurring when a <code>g</code>
|
||||
* switching function sign changes.</p>
|
||||
*
|
||||
* <p>Since events are only problem-dependent and are triggered by the
|
||||
* independent <i>time</i> variable and the state vector, they can
|
||||
* occur at virtually any time, unknown in advance. The integrators will
|
||||
* take care to avoid sign changes inside the steps, they will reduce
|
||||
* the step size when such an event is detected in order to put this
|
||||
* event exactly at the end of the current step. This guarantees that
|
||||
* step interpolation (which always has a one step scope) is relevant
|
||||
* even in presence of discontinuities. This is independent from the
|
||||
* stepsize control provided by integrators that monitor the local
|
||||
* error (this event handling feature is available for all integrators,
|
||||
* including fixed step ones).</p>
|
||||
*
|
||||
* @version $Revision$ $Date$
|
||||
* @since 2.1
|
||||
*/
|
||||
|
||||
public interface EventHandlerWithJacobians {
|
||||
|
||||
/** Stop indicator.
|
||||
* <p>This value should be used as the return value of the {@link
|
||||
* #eventOccurred eventOccurred} method when the integration should be
|
||||
* stopped after the event ending the current step.</p>
|
||||
*/
|
||||
int STOP = 0;
|
||||
|
||||
/** Reset state indicator.
|
||||
* <p>This value should be used as the return value of the {@link
|
||||
* #eventOccurred eventOccurred} method when the integration should
|
||||
* go on after the event ending the current step, with a new state
|
||||
* vector (which will be retrieved thanks to the {@link #resetState
|
||||
* resetState} method).</p>
|
||||
*/
|
||||
int RESET_STATE = 1;
|
||||
|
||||
/** Reset derivatives indicator.
|
||||
* <p>This value should be used as the return value of the {@link
|
||||
* #eventOccurred eventOccurred} method when the integration should
|
||||
* go on after the event ending the current step, with a new derivatives
|
||||
* vector (which will be retrieved thanks to the {@link
|
||||
* org.apache.commons.math.ode.FirstOrderDifferentialEquations#computeDerivatives}
|
||||
* method).</p>
|
||||
*/
|
||||
int RESET_DERIVATIVES = 2;
|
||||
|
||||
/** Continue indicator.
|
||||
* <p>This value should be used as the return value of the {@link
|
||||
* #eventOccurred eventOccurred} method when the integration should go
|
||||
* on after the event ending the current step.</p>
|
||||
*/
|
||||
int CONTINUE = 3;
|
||||
|
||||
/** Compute the value of the switching function.
|
||||
|
||||
* <p>The discrete events are generated when the sign of this
|
||||
* switching function changes. The integrator will take care to change
|
||||
* the stepsize in such a way these events occur exactly at step boundaries.
|
||||
* The switching function must be continuous in its roots neighborhood
|
||||
* (but not necessarily smooth), as the integrator will need to find its
|
||||
* roots to locate precisely the events.</p>
|
||||
|
||||
* @param t current value of the independent <i>time</i> variable
|
||||
* @param y array containing the current value of the state vector
|
||||
* @param dydy0 array containing the current value of the jacobian of
|
||||
* the state vector with respect to initial state
|
||||
* @param dydp array containing the current value of the jacobian of
|
||||
* the state vector with respect to parameters
|
||||
* @return value of the g switching function
|
||||
* @exception EventException if the switching function cannot be evaluated
|
||||
*/
|
||||
double g(double t, double[] y, double[][] dydy0, double[][] dydp)
|
||||
throws EventException;
|
||||
|
||||
/** Handle an event and choose what to do next.
|
||||
|
||||
* <p>This method is called when the integrator has accepted a step
|
||||
* ending exactly on a sign change of the function, just <em>before</em>
|
||||
* the step handler itself is called (see below for scheduling). It
|
||||
* allows the user to update his internal data to acknowledge the fact
|
||||
* the event has been handled (for example setting a flag in the {@link
|
||||
* org.apache.commons.math.ode.jacobians.ParameterizedODEWithJacobians
|
||||
* differential equations} to switch the derivatives computation in
|
||||
* case of discontinuity), or to direct the integrator to either stop
|
||||
* or continue integration, possibly with a reset state or derivatives.</p>
|
||||
|
||||
* <ul>
|
||||
* <li>if {@link #STOP} is returned, the step handler will be called
|
||||
* with the <code>isLast</code> flag of the {@link
|
||||
* org.apache.commons.math.ode.jacobians.StepHandlerWithJacobians#handleStep(
|
||||
* StepInterpolatorWithJacobians, boolean) handleStep} method set to true and
|
||||
* the integration will be stopped,</li>
|
||||
* <li>if {@link #RESET_STATE} is returned, the {@link #resetState
|
||||
* resetState} method will be called once the step handler has
|
||||
* finished its task, and the integrator will also recompute the
|
||||
* derivatives,</li>
|
||||
* <li>if {@link #RESET_DERIVATIVES} is returned, the integrator
|
||||
* will recompute the derivatives,
|
||||
* <li>if {@link #CONTINUE} is returned, no specific action will
|
||||
* be taken (apart from having called this method) and integration
|
||||
* will continue.</li>
|
||||
* </ul>
|
||||
|
||||
* <p>The scheduling between this method and the {@link
|
||||
* org.apache.commons.math.ode.jacobians.StepHandlerWithJacobians
|
||||
* StepHandlerWithJacobians} method {@link
|
||||
* org.apache.commons.math.ode.jacobians.StepHandlerWithJacobians#handleStep(
|
||||
* StepInterpolatorWithJacobians, boolean) handleStep(interpolator, isLast)}
|
||||
* is to call this method first and <code>handleStep</code> afterwards. This
|
||||
* scheduling allows the integrator to pass <code>true</code> as the
|
||||
* <code>isLast</code> parameter to the step handler to make it aware the step
|
||||
* will be the last one if this method returns {@link #STOP}. As the
|
||||
* interpolator may be used to navigate back throughout the last step (as {@link
|
||||
* org.apache.commons.math.ode.sampling.StepNormalizer StepNormalizer}
|
||||
* does for example), user code called by this method and user
|
||||
* code called by step handlers may experience apparently out of order values
|
||||
* of the independent time variable. As an example, if the same user object
|
||||
* implements both this {@link EventHandlerWithJacobians EventHandler} interface and the
|
||||
* {@link org.apache.commons.math.ode.sampling.FixedStepHandler FixedStepHandler}
|
||||
* interface, a <em>forward</em> integration may call its
|
||||
* <code>eventOccurred</code> method with t = 10 first and call its
|
||||
* <code>handleStep</code> method with t = 9 afterwards. Such out of order
|
||||
* calls are limited to the size of the integration step for {@link
|
||||
* org.apache.commons.math.ode.sampling.StepHandler variable step handlers} and
|
||||
* to the size of the fixed step for {@link
|
||||
* org.apache.commons.math.ode.sampling.FixedStepHandler fixed step handlers}.</p>
|
||||
|
||||
* @param t current value of the independent <i>time</i> variable
|
||||
* @param y array containing the current value of the state vector
|
||||
* @param dydy0 array containing the current value of the jacobian of
|
||||
* the state vector with respect to initial state
|
||||
* @param dydp array containing the current value of the jacobian of
|
||||
* the state vector with respect to parameters
|
||||
* @param increasing if true, the value of the switching function increases
|
||||
* when times increases around event (note that increase is measured with respect
|
||||
* to physical time, not with respect to integration which may go backward in time)
|
||||
* @return indication of what the integrator should do next, this
|
||||
* value must be one of {@link #STOP}, {@link #RESET_STATE},
|
||||
* {@link #RESET_DERIVATIVES} or {@link #CONTINUE}
|
||||
* @exception EventException if the event occurrence triggers an error
|
||||
*/
|
||||
int eventOccurred(double t, double[] y, double[][] dydy0, double[][] dydp,
|
||||
boolean increasing) throws EventException;
|
||||
|
||||
/** Reset the state prior to continue the integration.
|
||||
|
||||
* <p>This method is called after the step handler has returned and
|
||||
* before the next step is started, but only when {@link
|
||||
* #eventOccurred} has itself returned the {@link #RESET_STATE}
|
||||
* indicator. It allows the user to reset the state vector for the
|
||||
* next step, without perturbing the step handler of the finishing
|
||||
* step. If the {@link #eventOccurred} never returns the {@link
|
||||
* #RESET_STATE} indicator, this function will never be called, and it is
|
||||
* safe to leave its body empty.</p>
|
||||
|
||||
* @param t current value of the independent <i>time</i> variable
|
||||
* @param y array containing the current value of the state vector
|
||||
* the new state should be put in the same array
|
||||
* @param dydy0 array containing the current value of the jacobian of
|
||||
* the state vector with respect to initial state, the new jacobian
|
||||
* should be put in the same array
|
||||
* @param dydp array containing the current value of the jacobian of
|
||||
* the state vector with respect to parameters, the new jacobian
|
||||
* should be put in the same array
|
||||
* @exception EventException if the state cannot be reseted
|
||||
*/
|
||||
void resetState(double t, double[] y, double[][] dydy0, double[][] dydp)
|
||||
throws EventException;
|
||||
|
||||
}
|
||||
@ -30,6 +30,8 @@ import org.apache.commons.math.ode.DerivativeException;
|
||||
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
|
||||
import org.apache.commons.math.ode.FirstOrderIntegrator;
|
||||
import org.apache.commons.math.ode.IntegratorException;
|
||||
import org.apache.commons.math.ode.events.EventException;
|
||||
import org.apache.commons.math.ode.events.EventHandler;
|
||||
import org.apache.commons.math.ode.sampling.StepHandler;
|
||||
import org.apache.commons.math.ode.sampling.StepInterpolator;
|
||||
|
||||
@ -129,23 +131,50 @@ public class FirstOrderIntegratorWithJacobians {
|
||||
* @see #getStepHandlers()
|
||||
*/
|
||||
public void clearStepHandlers() {
|
||||
integrator.clearStepHandlers();
|
||||
}
|
||||
|
||||
// preserve the handlers we did not add ourselves
|
||||
final Collection<StepHandler> otherHandlers = new ArrayList<StepHandler>();
|
||||
for (final StepHandler handler : integrator.getStepHandlers()) {
|
||||
if (!(handler instanceof StepHandlerWrapper)) {
|
||||
otherHandlers.add(handler);
|
||||
/** Add an event handler to the integrator.
|
||||
* @param handler event handler
|
||||
* @param maxCheckInterval maximal time interval between switching
|
||||
* function checks (this interval prevents missing sign changes in
|
||||
* case the integration steps becomes very large)
|
||||
* @param convergence convergence threshold in the event time search
|
||||
* @param maxIterationCount upper limit of the iteration count in
|
||||
* the event time search
|
||||
* @see #getEventHandlers()
|
||||
* @see #clearEventHandlers()
|
||||
*/
|
||||
public void addEventHandler(EventHandlerWithJacobians handler,
|
||||
double maxCheckInterval,
|
||||
double convergence,
|
||||
int maxIterationCount) {
|
||||
integrator.addEventHandler(new EventHandlerWrapper(handler),
|
||||
maxCheckInterval, convergence, maxIterationCount);
|
||||
}
|
||||
|
||||
/** Get all the event handlers that have been added to the integrator.
|
||||
* @return an unmodifiable collection of the added events handlers
|
||||
* @see #addEventHandler(EventHandlerWithJacobians, double, double, int)
|
||||
* @see #clearEventHandlers()
|
||||
*/
|
||||
public Collection<EventHandlerWithJacobians> getEventHandlers() {
|
||||
final Collection<EventHandlerWithJacobians> handlers =
|
||||
new ArrayList<EventHandlerWithJacobians>();
|
||||
for (final EventHandler handler : integrator.getEventHandlers()) {
|
||||
if (handler instanceof EventHandlerWrapper) {
|
||||
handlers.add(((EventHandlerWrapper) handler).getHandler());
|
||||
}
|
||||
}
|
||||
return handlers;
|
||||
}
|
||||
|
||||
// clear all handlers
|
||||
integrator.clearStepHandlers();
|
||||
|
||||
// put back the preserved handlers
|
||||
for (final StepHandler handler : otherHandlers) {
|
||||
integrator.addStepHandler(handler);
|
||||
}
|
||||
|
||||
/** Remove all the event handlers that have been added to the integrator.
|
||||
* @see #addEventHandler(EventHandlerWithJacobians, double, double, int)
|
||||
* @see #getEventHandlers()
|
||||
*/
|
||||
public void clearEventHandlers() {
|
||||
integrator.clearEventHandlers();
|
||||
}
|
||||
|
||||
/** Integrate the differential equations and the variational equations up to the given time.
|
||||
@ -217,18 +246,39 @@ public class FirstOrderIntegratorWithJacobians {
|
||||
final double stopTime = integrator.integrate(new MappingWrapper(), t0, z, t, z);
|
||||
|
||||
// dispatch the final compound state into the state and partial derivatives arrays
|
||||
System.arraycopy(z, 0, y, 0, n);
|
||||
for (int i = 0; i < n; ++i) {
|
||||
System.arraycopy(z, n * (i + 1), dYdY0[i], 0, n);
|
||||
}
|
||||
for (int i = 0; i < n; ++i) {
|
||||
System.arraycopy(z, n * (n + 1) + i * k, dYdP[i], 0, k);
|
||||
}
|
||||
dispatchCompoundState(z, y, dYdY0, dYdP);
|
||||
|
||||
return stopTime;
|
||||
|
||||
}
|
||||
|
||||
/** Dispatch a compound state array into state and jacobians arrays.
|
||||
* @param z compound state
|
||||
* @param y raw state array to fill
|
||||
* @param dydy0 jacobian array to fill
|
||||
* @param dydp jacobian array to fill
|
||||
*/
|
||||
private static void dispatchCompoundState(final double[] z, final double[] y,
|
||||
final double[][] dydy0, final double[][] dydp) {
|
||||
|
||||
final int n = y.length;
|
||||
final int k = dydp[0].length;
|
||||
|
||||
// state
|
||||
System.arraycopy(z, 0, y, 0, n);
|
||||
|
||||
// jacobian with respect to initial state
|
||||
for (int i = 0; i < n; ++i) {
|
||||
System.arraycopy(z, n * (i + 1), dydy0[i], 0, n);
|
||||
}
|
||||
|
||||
// jacobian with respect to parameters
|
||||
for (int i = 0; i < n; ++i) {
|
||||
System.arraycopy(z, n * (n + 1) + i * k, dydp[i], 0, k);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/** Get the current value of the step start time t<sub>i</sub>.
|
||||
* <p>This method can be called during integration (typically by
|
||||
* the object implementing the {@link FirstOrderDifferentialEquations
|
||||
@ -766,4 +816,62 @@ public class FirstOrderIntegratorWithJacobians {
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/** Wrapper for event handlers. */
|
||||
private class EventHandlerWrapper implements EventHandler {
|
||||
|
||||
/** Underlying event handler with jacobians. */
|
||||
private final EventHandlerWithJacobians handler;
|
||||
|
||||
/** State array. */
|
||||
private double[] y;
|
||||
|
||||
/** Jacobian with respect to initial state dy/dy0. */
|
||||
private double[][] dydy0;
|
||||
|
||||
/** Jacobian with respect to parameters dy/dp. */
|
||||
private double[][] dydp;
|
||||
|
||||
/** Simple constructor.
|
||||
* @param handler underlying event handler with jacobians
|
||||
*/
|
||||
public EventHandlerWrapper(final EventHandlerWithJacobians handler) {
|
||||
this.handler = handler;
|
||||
final int n = ode.getDimension();
|
||||
final int k = ode.getParametersDimension();
|
||||
y = new double[n];
|
||||
dydy0 = new double[n][n];
|
||||
dydp = new double[n][k];
|
||||
}
|
||||
|
||||
/** Get the underlying event handler with jacobians.
|
||||
* @return underlying event handler with jacobians
|
||||
*/
|
||||
public EventHandlerWithJacobians getHandler() {
|
||||
return handler;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public int eventOccurred(double t, double[] z, boolean increasing)
|
||||
throws EventException {
|
||||
dispatchCompoundState(z, y, dydy0, dydp);
|
||||
return handler.eventOccurred(t, y, dydy0, dydp, increasing);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double g(double t, double[] z)
|
||||
throws EventException {
|
||||
dispatchCompoundState(z, y, dydy0, dydp);
|
||||
return handler.g(t, y, dydy0, dydp);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public void resetState(double t, double[] z)
|
||||
throws EventException {
|
||||
dispatchCompoundState(z, y, dydy0, dydp);
|
||||
handler.resetState(t, y, dydy0, dydp);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@ -242,6 +242,39 @@ public class FirstOrderIntegratorWithJacobiansTest {
|
||||
extInt.integrate(0, y, circle.exactDyDp(0), t, y, dydy0, dydp);
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testEventHandler() throws IntegratorException, DerivativeException {
|
||||
FirstOrderIntegrator integ =
|
||||
new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10);
|
||||
double[] y = new double[] { 0.0, 1.0 };
|
||||
final Circle circle = new Circle(y, 1.0, 1.0, 0.1);
|
||||
double[][] dydy0 = new double[2][2];
|
||||
double[][] dydp = new double[2][3];
|
||||
double t = 18 * Math.PI;
|
||||
final FirstOrderIntegratorWithJacobians extInt =
|
||||
new FirstOrderIntegratorWithJacobians(integ, circle);
|
||||
extInt.addEventHandler(new EventHandlerWithJacobians() {
|
||||
|
||||
public int eventOccurred(double t, double[] y, double[][] dydy0,
|
||||
double[][] dydp, boolean increasing) {
|
||||
Assert.assertEquals(0.1, y[1], 1.0e-11);
|
||||
Assert.assertTrue(!increasing);
|
||||
return STOP;
|
||||
}
|
||||
|
||||
public double g(double t, double[] y, double[][] dydy0,
|
||||
double[][] dydp) {
|
||||
return y[1] - 0.1;
|
||||
}
|
||||
|
||||
public void resetState(double t, double[] y, double[][] dydy0,
|
||||
double[][] dydp) {
|
||||
}
|
||||
}, 10.0, 1.0e-10, 1000);
|
||||
double stopTime = extInt.integrate(0, y, circle.exactDyDp(0), t, y, dydy0, dydp);
|
||||
Assert.assertTrue(stopTime < 5.0 * Math.PI);
|
||||
}
|
||||
|
||||
private static class Brusselator implements ParameterizedODEWithJacobians {
|
||||
|
||||
private double b;
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user