Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

removed the ode.jacobians package 

Jira: MATH-380

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1037343 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2010-11-20 22:02:55 +00:00
parent 3859aa17cf
commit 22d63af83d
9 changed files with 8 additions and 2172 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

224
src/main/java/org/apache/commons/math/ode/jacobians/EventHandlerWithJacobians.java
View File

@ -1,224 +0,0 @@
/*
* 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>
*
* <p>Note that is is possible to register a {@link
* org.apache.commons.math.ode.events.EventHandler classical event handler}
* in the low level integrator used to build a {@link FirstOrderIntegratorWithJacobians}
* rather than implementing this class. The event handlers registered at low level
* will see the big compound state whether the event handlers defined by this interface
* see the original state, and its jacobians in separate arrays.</p>
*
* <p>The compound state is guaranteed to contain the original state in the first
* elements, followed by the jacobian with respect to initial state (in row order),
* followed by the jacobian with respect to parameters (in row order). If for example
* the original state dimension is 6 and there are 3 parameters, the compound state will
* be a 60 elements array. The first 6 elements will be the original state, the next 36
* elements will be the jacobian with respect to initial state, and the remaining 18 elements
* will be the jacobian with respect to parameters.</p>
*
* <p>Dealing with low level event handlers is cumbersome if one really needs the jacobians
* in these methods, but it also prevents many data being copied back and forth between
* state and jacobians on one side and compound state on the other side. So for performance
* reasons, it is recommended to use this interface <em>only</em> if jacobians are really
* needed and to use lower level handlers if only state is needed.</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.ODEWithJacobians
* 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;
}

896
src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java
View File

@ -1,896 +0,0 @@
/*
* 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 java.io.IOException;
import java.io.ObjectInput;
import java.io.ObjectOutput;
import java.lang.reflect.Array;
import java.util.ArrayList;
import java.util.Collection;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.MaxEvaluationsExceededException;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.ode.ExtendedFirstOrderDifferentialEquations;
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;
/** This class enhances a first order integrator for differential equations to
* compute also partial derivatives of the solution with respect to initial state
* and parameters.
* <p>In order to compute both the state and its derivatives, the ODE problem
* is extended with jacobians of the raw ODE and the variational equations are
* added to form a new compound problem of higher dimension. If the original ODE
* problem has dimension n and there are p parameters, the compound problem will
* have dimension n &times; (1 + n + p).</p>
* @see ParameterizedODE
* @see ODEWithJacobians
* @version $Revision$ $Date$
* @since 2.1
*/
public class FirstOrderIntegratorWithJacobians {
/** Underlying integrator for compound problem. */
private final FirstOrderIntegrator integrator;
/** Raw equations to integrate. */
private final ODEWithJacobians ode;
/** Maximal number of evaluations allowed. */
private int maxEvaluations;
/** Number of evaluations already performed. */
private int evaluations;
/** Build an enhanced integrator using internal differentiation to compute jacobians.
* @param integrator underlying integrator to solve the compound problem
* @param ode original problem (f in the equation y' = f(t, y))
* @param p parameters array (may be null if {@link
* ParameterizedODE#getParametersDimension()
* getParametersDimension()} from original problem is zero)
* @param hY step sizes to use for computing the jacobian df/dy, must have the
* same dimension as the original problem
* @param hP step sizes to use for computing the jacobian df/dp, must have the
* same dimension as the original problem parameters dimension
* @see #FirstOrderIntegratorWithJacobians(FirstOrderIntegrator,
* ODEWithJacobians)
*/
public FirstOrderIntegratorWithJacobians(final FirstOrderIntegrator integrator,
final ParameterizedODE ode,
final double[] p, final double[] hY, final double[] hP) {
checkDimension(ode.getDimension(), hY);
checkDimension(ode.getParametersDimension(), p);
checkDimension(ode.getParametersDimension(), hP);
this.integrator = integrator;
this.ode = new FiniteDifferencesWrapper(ode, p, hY, hP);
setMaxEvaluations(-1);
}
/** Build an enhanced integrator using ODE builtin jacobian computation features.
* @param integrator underlying integrator to solve the compound problem
* @param ode original problem, which can compute the jacobians by itself
* @see #FirstOrderIntegratorWithJacobians(FirstOrderIntegrator,
* ParameterizedODE, double[], double[], double[])
*/
public FirstOrderIntegratorWithJacobians(final FirstOrderIntegrator integrator,
final ODEWithJacobians ode) {
this.integrator = integrator;
this.ode = ode;
setMaxEvaluations(-1);
}
/** Add a step handler to this integrator.
* <p>The handler will be called by the integrator for each accepted
* step.</p>
* @param handler handler for the accepted steps
* @see #getStepHandlers()
* @see #clearStepHandlers()
*/
public void addStepHandler(StepHandlerWithJacobians handler) {
final int n = ode.getDimension();
final int k = ode.getParametersDimension();
integrator.addStepHandler(new StepHandlerWrapper(handler, n, k));
}
/** Get all the step handlers that have been added to the integrator.
* @return an unmodifiable collection of the added events handlers
* @see #addStepHandler(StepHandlerWithJacobians)
* @see #clearStepHandlers()
*/
public Collection<StepHandlerWithJacobians> getStepHandlers() {
final Collection<StepHandlerWithJacobians> handlers =
new ArrayList<StepHandlerWithJacobians>();
for (final StepHandler handler : integrator.getStepHandlers()) {
if (handler instanceof StepHandlerWrapper) {
handlers.add(((StepHandlerWrapper) handler).getHandler());
}
}
return handlers;
}
/** Remove all the step handlers that have been added to the integrator.
* @see #addStepHandler(StepHandlerWithJacobians)
* @see #getStepHandlers()
*/
public void clearStepHandlers() {
integrator.clearStepHandlers();
}
/** 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) {
final int n = ode.getDimension();
final int k = ode.getParametersDimension();
integrator.addEventHandler(new EventHandlerWrapper(handler, n, k),
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;
}
/** 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.
* <p>This method solves an Initial Value Problem (IVP) and also computes the derivatives
* of the solution with respect to initial state and parameters. This can be used as
* a basis to solve Boundary Value Problems (BVP).</p>
* <p>Since this method stores some internal state variables made
* available in its public interface during integration ({@link
* #getCurrentSignedStepsize()}), it is <em>not</em> thread-safe.</p>
* @param t0 initial time
* @param y0 initial value of the state vector at t0
* @param dY0dP initial value of the state vector derivative with respect to the
* parameters at t0
* @param t target time for the integration
* (can be set to a value smaller than <code>t0</code> for backward integration)
* @param y placeholder where to put the state vector at each successful
* step (and hence at the end of integration), can be the same object as y0
* @param dYdY0 placeholder where to put the state vector derivative with respect
* to the initial state (dy[i]/dy0[j] is in element array dYdY0[i][j]) at each successful
* step (and hence at the end of integration)
* @param dYdP placeholder where to put the state vector derivative with respect
* to the parameters (dy[i]/dp[j] is in element array dYdP[i][j]) at each successful
* step (and hence at the end of integration)
* @return stop time, will be the same as target time if integration reached its
* target, but may be different if some event handler stops it at some point.
* @throws IntegratorException if the integrator cannot perform integration
* @throws MathUserException this exception is propagated to the caller if
* the underlying user function triggers one
*/
public double integrate(final double t0, final double[] y0, final double[][] dY0dP,
final double t, final double[] y,
final double[][] dYdY0, final double[][] dYdP)
throws MathUserException, IntegratorException {
final int n = ode.getDimension();
final int k = ode.getParametersDimension();
checkDimension(n, y0);
checkDimension(n, y);
checkDimension(n, dYdY0);
checkDimension(n, dYdY0[0]);
if (k != 0) {
checkDimension(n, dY0dP);
checkDimension(k, dY0dP[0]);
checkDimension(n, dYdP);
checkDimension(k, dYdP[0]);
}
// set up initial state, including partial derivatives
// the compound state z contains the raw state y and its derivatives
// with respect to initial state y0 and to parameters p
// y[i] is stored in z[i]
// dy[i]/dy0[j] is stored in z[n + i * n + j]
// dy[i]/dp[j] is stored in z[n * (n + 1) + i * k + j]
final double[] z = new double[n * (1 + n + k)];
System.arraycopy(y0, 0, z, 0, n);
for (int i = 0; i < n; ++i) {
// set diagonal element of dy/dy0 to 1.0 at t = t0
z[i * (1 + n) + n] = 1.0;
// set initial derivatives with respect to parameters
System.arraycopy(dY0dP[i], 0, z, n * (n + 1) + i * k, k);
}
// integrate the compound state variational equations
evaluations = 0;
final double stopTime = integrator.integrate(new MappingWrapper(), t0, z, t, z);
// dispatch the final compound state into the state and partial derivatives arrays
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 org.apache.commons.math.ode.FirstOrderDifferentialEquations
* differential equations} problem) if the value of the current step that
* is attempted is needed.</p>
* <p>The result is undefined if the method is called outside of
* calls to <code>integrate</code>.</p>
* @return current value of the step start time t<sub>i</sub>
*/
public double getCurrentStepStart() {
return integrator.getCurrentStepStart();
}
/** Get the current signed value of the integration stepsize.
* <p>This method can be called during integration (typically by
* the object implementing the {@link org.apache.commons.math.ode.FirstOrderDifferentialEquations
* differential equations} problem) if the signed value of the current stepsize
* that is tried is needed.</p>
* <p>The result is undefined if the method is called outside of
* calls to <code>integrate</code>.</p>
* @return current signed value of the stepsize
*/
public double getCurrentSignedStepsize() {
return integrator.getCurrentSignedStepsize();
}
/** Set the maximal number of differential equations function evaluations.
* <p>The purpose of this method is to avoid infinite loops which can occur
* for example when stringent error constraints are set or when lots of
* discrete events are triggered, thus leading to many rejected steps.</p>
* @param maxEvaluations maximal number of function evaluations (negative
* values are silently converted to maximal integer value, thus representing
* almost unlimited evaluations)
*/
public void setMaxEvaluations(int maxEvaluations) {
this.maxEvaluations = (maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations;
}
/** Get the maximal number of functions evaluations.
* @return maximal number of functions evaluations
*/
public int getMaxEvaluations() {
return maxEvaluations;
}
/** Get the number of evaluations of the differential equations function.
* <p>
* The number of evaluations corresponds to the last call to the
* <code>integrate</code> method. It is 0 if the method has not been called yet.
* </p>
* @return number of evaluations of the differential equations function
*/
public int getEvaluations() {
return evaluations;
}
/** Check array dimensions.
* @param expected expected dimension
* @param array (may be null if expected is 0)
* @throws IllegalArgumentException if the array dimension does not match the expected one
*/
private void checkDimension(final int expected, final Object array)
throws IllegalArgumentException {
int arrayDimension = (array == null) ? 0 : Array.getLength(array);
if (arrayDimension != expected) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, arrayDimension, expected);
}
}
/** Wrapper class used to map state and jacobians into compound state. */
private class MappingWrapper implements ExtendedFirstOrderDifferentialEquations {
/** Current state. */
private final double[] y;
/** Time derivative of the current state. */
private final double[] yDot;
/** Derivatives of yDot with respect to state. */
private final double[][] dFdY;
/** Derivatives of yDot with respect to parameters. */
private final double[][] dFdP;
/** Simple constructor.
*/
public MappingWrapper() {
final int n = ode.getDimension();
final int k = ode.getParametersDimension();
y = new double[n];
yDot = new double[n];
dFdY = new double[n][n];
dFdP = new double[n][k];
}
/** {@inheritDoc} */
public int getDimension() {
final int n = y.length;
final int k = dFdP[0].length;
return n * (1 + n + k);
}
/** {@inheritDoc} */
public int getMainSetDimension() {
return ode.getDimension();
}
/** {@inheritDoc} */
public void computeDerivatives(final double t, final double[] z, final double[] zDot)
throws MathUserException {
final int n = y.length;
final int k = dFdP[0].length;
// compute raw ODE and its jacobians: dy/dt, d[dy/dt]/dy0 and d[dy/dt]/dp
System.arraycopy(z, 0, y, 0, n);
if (++evaluations > maxEvaluations) {
throw new MathUserException(new MaxEvaluationsExceededException(maxEvaluations));
}
ode.computeDerivatives(t, y, yDot);
ode.computeJacobians(t, y, yDot, dFdY, dFdP);
// state part of the compound equations
System.arraycopy(yDot, 0, zDot, 0, n);
// variational equations: from d[dy/dt]/dy0 to d[dy/dy0]/dt
for (int i = 0; i < n; ++i) {
final double[] dFdYi = dFdY[i];
for (int j = 0; j < n; ++j) {
double s = 0;
final int startIndex = n + j;
int zIndex = startIndex;
for (int l = 0; l < n; ++l) {
s += dFdYi[l] * z[zIndex];
zIndex += n;
}
zDot[startIndex + i * n] = s;
}
}
// variational equations: from d[dy/dt]/dy0 and d[dy/dt]/dp to d[dy/dp]/dt
for (int i = 0; i < n; ++i) {
final double[] dFdYi = dFdY[i];
final double[] dFdPi = dFdP[i];
for (int j = 0; j < k; ++j) {
double s = dFdPi[j];
final int startIndex = n * (n + 1) + j;
int zIndex = startIndex;
for (int l = 0; l < n; ++l) {
s += dFdYi[l] * z[zIndex];
zIndex += k;
}
zDot[startIndex + i * k] = s;
}
}
}
}
/** Wrapper class to compute jacobians by finite differences for ODE which do not compute them themselves. */
private class FiniteDifferencesWrapper implements ODEWithJacobians {
/** Raw ODE without jacobians computation. */
private final ParameterizedODE ode;
/** Parameters array (may be null if parameters dimension from original problem is zero) */
private final double[] p;
/** Step sizes to use for computing the jacobian df/dy. */
private final double[] hY;
/** Step sizes to use for computing the jacobian df/dp. */
private final double[] hP;
/** Temporary array for state derivatives used to compute jacobians. */
private final double[] tmpDot;
/** Simple constructor.
* @param ode original ODE problem, without jacobians computations
* @param p parameters array (may be null if parameters dimension from original problem is zero)
* @param hY step sizes to use for computing the jacobian df/dy
* @param hP step sizes to use for computing the jacobian df/dp
*/
public FiniteDifferencesWrapper(final ParameterizedODE ode,
final double[] p, final double[] hY, final double[] hP) {
this.ode = ode;
this.p = p.clone();
this.hY = hY.clone();
this.hP = hP.clone();
tmpDot = new double[ode.getDimension()];
}
/** {@inheritDoc} */
public int getDimension() {
return ode.getDimension();
}
/** {@inheritDoc} */
public void computeDerivatives(double t, double[] y, double[] yDot) throws MathUserException {
// this call to computeDerivatives has already been counted,
// we must not increment the counter again
ode.computeDerivatives(t, y, yDot);
}
/** {@inheritDoc} */
public int getParametersDimension() {
return ode.getParametersDimension();
}
/** {@inheritDoc} */
public void computeJacobians(double t, double[] y, double[] yDot,
double[][] dFdY, double[][] dFdP)
throws MathUserException {
final int n = hY.length;
final int k = hP.length;
evaluations += n + k;
if (evaluations > maxEvaluations) {
throw new MathUserException(new MaxEvaluationsExceededException(maxEvaluations));
}
// compute df/dy where f is the ODE and y is the state array
for (int j = 0; j < n; ++j) {
final double savedYj = y[j];
y[j] += hY[j];
ode.computeDerivatives(t, y, tmpDot);
for (int i = 0; i < n; ++i) {
dFdY[i][j] = (tmpDot[i] - yDot[i]) / hY[j];
}
y[j] = savedYj;
}
// compute df/dp where f is the ODE and p is the parameters array
for (int j = 0; j < k; ++j) {
ode.setParameter(j, p[j] + hP[j]);
ode.computeDerivatives(t, y, tmpDot);
for (int i = 0; i < n; ++i) {
dFdP[i][j] = (tmpDot[i] - yDot[i]) / hP[j];
}
ode.setParameter(j, p[j]);
}
}
}
/** Wrapper for step handlers. */
private static class StepHandlerWrapper implements StepHandler {
/** Underlying step handler with jacobians. */
private final StepHandlerWithJacobians handler;
/** Dimension of the original ODE. */
private final int n;
/** Number of parameters. */
private final int k;
/** Simple constructor.
* @param handler underlying step handler with jacobians
* @param n dimension of the original ODE
* @param k number of parameters
*/
public StepHandlerWrapper(final StepHandlerWithJacobians handler,
final int n, final int k) {
this.handler = handler;
this.n = n;
this.k = k;
}
/** Get the underlying step handler with jacobians.
* @return underlying step handler with jacobians
*/
public StepHandlerWithJacobians getHandler() {
return handler;
}
/** {@inheritDoc} */
public void handleStep(StepInterpolator interpolator, boolean isLast)
throws MathUserException {
handler.handleStep(new StepInterpolatorWrapper(interpolator, n, k), isLast);
}
/** {@inheritDoc} */
public boolean requiresDenseOutput() {
return handler.requiresDenseOutput();
}
/** {@inheritDoc} */
public void reset() {
handler.reset();
}
}
/** Wrapper for step interpolators. */
private static class StepInterpolatorWrapper
implements StepInterpolatorWithJacobians {
/** Wrapped interpolator. */
private StepInterpolator interpolator;
/** 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;
/** Time derivative of the state array. */
private double[] yDot;
/** Time derivative of the sacobian with respect to initial state dy/dy0. */
private double[][] dydy0Dot;
/** Time derivative of the jacobian with respect to parameters dy/dp. */
private double[][] dydpDot;
/** Simple constructor.
* <p>This constructor is used only for externalization. It does nothing.</p>
*/
@SuppressWarnings("unused")
public StepInterpolatorWrapper() {
}
/** Simple constructor.
* @param interpolator wrapped interpolator
* @param n dimension of the original ODE
* @param k number of parameters
*/
public StepInterpolatorWrapper(final StepInterpolator interpolator,
final int n, final int k) {
this.interpolator = interpolator;
y = new double[n];
dydy0 = new double[n][n];
dydp = new double[n][k];
yDot = new double[n];
dydy0Dot = new double[n][n];
dydpDot = new double[n][k];
}
/** {@inheritDoc} */
public void setInterpolatedTime(double time) {
interpolator.setInterpolatedTime(time);
}
/** {@inheritDoc} */
public boolean isForward() {
return interpolator.isForward();
}
/** {@inheritDoc} */
public double getPreviousTime() {
return interpolator.getPreviousTime();
}
/** {@inheritDoc} */
public double getInterpolatedTime() {
return interpolator.getInterpolatedTime();
}
/** {@inheritDoc} */
public double[] getInterpolatedY() throws MathUserException {
double[] extendedState = interpolator.getInterpolatedState();
System.arraycopy(extendedState, 0, y, 0, y.length);
return y;
}
/** {@inheritDoc} */
public double[][] getInterpolatedDyDy0() throws MathUserException {
double[] extendedState = interpolator.getInterpolatedState();
final int n = y.length;
int start = n;
for (int i = 0; i < n; ++i) {
System.arraycopy(extendedState, start, dydy0[i], 0, n);
start += n;
}
return dydy0;
}
/** {@inheritDoc} */
public double[][] getInterpolatedDyDp() throws MathUserException {
double[] extendedState = interpolator.getInterpolatedState();
final int n = y.length;
final int k = dydp[0].length;
int start = n * (n + 1);
for (int i = 0; i < n; ++i) {
System.arraycopy(extendedState, start, dydp[i], 0, k);
start += k;
}
return dydp;
}
/** {@inheritDoc} */
public double[] getInterpolatedYDot() throws MathUserException {
double[] extendedDerivatives = interpolator.getInterpolatedDerivatives();
System.arraycopy(extendedDerivatives, 0, yDot, 0, yDot.length);
return yDot;
}
/** {@inheritDoc} */
public double[][] getInterpolatedDyDy0Dot() throws MathUserException {
double[] extendedDerivatives = interpolator.getInterpolatedDerivatives();
final int n = y.length;
int start = n;
for (int i = 0; i < n; ++i) {
System.arraycopy(extendedDerivatives, start, dydy0Dot[i], 0, n);
start += n;
}
return dydy0Dot;
}
/** {@inheritDoc} */
public double[][] getInterpolatedDyDpDot() throws MathUserException {
double[] extendedDerivatives = interpolator.getInterpolatedDerivatives();
final int n = y.length;
final int k = dydpDot[0].length;
int start = n * (n + 1);
for (int i = 0; i < n; ++i) {
System.arraycopy(extendedDerivatives, start, dydpDot[i], 0, k);
start += k;
}
return dydpDot;
}
/** {@inheritDoc} */
public double getCurrentTime() {
return interpolator.getCurrentTime();
}
/** {@inheritDoc} */
public StepInterpolatorWithJacobians copy() throws MathUserException {
final int n = y.length;
final int k = dydp[0].length;
StepInterpolatorWrapper copied =
new StepInterpolatorWrapper(interpolator.copy(), n, k);
copyArray(y, copied.y);
copyArray(dydy0, copied.dydy0);
copyArray(dydp, copied.dydp);
copyArray(yDot, copied.yDot);
copyArray(dydy0Dot, copied.dydy0Dot);
copyArray(dydpDot, copied.dydpDot);
return copied;
}
/** {@inheritDoc} */
public void writeExternal(ObjectOutput out) throws IOException {
out.writeObject(interpolator);
out.writeInt(y.length);
out.writeInt(dydp[0].length);
writeArray(out, y);
writeArray(out, dydy0);
writeArray(out, dydp);
writeArray(out, yDot);
writeArray(out, dydy0Dot);
writeArray(out, dydpDot);
}
/** {@inheritDoc} */
public void readExternal(ObjectInput in) throws IOException, ClassNotFoundException {
interpolator = (StepInterpolator) in.readObject();
final int n = in.readInt();
final int k = in.readInt();
y = new double[n];
dydy0 = new double[n][n];
dydp = new double[n][k];
yDot = new double[n];
dydy0Dot = new double[n][n];
dydpDot = new double[n][k];
readArray(in, y);
readArray(in, dydy0);
readArray(in, dydp);
readArray(in, yDot);
readArray(in, dydy0Dot);
readArray(in, dydpDot);
}
/** Copy an array.
* @param src source array
* @param dest destination array
*/
private static void copyArray(final double[] src, final double[] dest) {
System.arraycopy(src, 0, dest, 0, src.length);
}
/** Copy an array.
* @param src source array
* @param dest destination array
*/
private static void copyArray(final double[][] src, final double[][] dest) {
for (int i = 0; i < src.length; ++i) {
copyArray(src[i], dest[i]);
}
}
/** Write an array.
* @param out output stream
* @param array array to write
* @exception IOException if array cannot be read
*/
private static void writeArray(final ObjectOutput out, final double[] array)
throws IOException {
for (int i = 0; i < array.length; ++i) {
out.writeDouble(array[i]);
}
}
/** Write an array.
* @param out output stream
* @param array array to write
* @exception IOException if array cannot be read
*/
private static void writeArray(final ObjectOutput out, final double[][] array)
throws IOException {
for (int i = 0; i < array.length; ++i) {
writeArray(out, array[i]);
}
}
/** Read an array.
* @param in input stream
* @param array array to read
* @exception IOException if array cannot be read
*/
private static void readArray(final ObjectInput in, final double[] array)
throws IOException {
for (int i = 0; i < array.length; ++i) {
array[i] = in.readDouble();
}
}
/** Read an array.
* @param in input stream
* @param array array to read
* @exception IOException if array cannot be read
*/
private static void readArray(final ObjectInput in, final double[][] array)
throws IOException {
for (int i = 0; i < array.length; ++i) {
readArray(in, array[i]);
}
}
}
/** Wrapper for event handlers. */
private static 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
* @param n dimension of the original ODE
* @param k number of parameters
*/
public EventHandlerWrapper(final EventHandlerWithJacobians handler,
final int n, final int k) {
this.handler = handler;
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);
}
}
}

52
src/main/java/org/apache/commons/math/ode/jacobians/ODEWithJacobians.java
View File

@ -1,52 +0,0 @@
/*
* 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.exception.MathUserException;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
/** This interface represents {@link ParameterizedODE
* first order differential equations} with parameters and partial derivatives.
*
* @see FirstOrderIntegratorWithJacobians
*
* @version $Revision$ $Date$
* @since 2.1
*/
public interface ODEWithJacobians extends FirstOrderDifferentialEquations {
/** Get the number of parameters.
* @return number of parameters
*/
int getParametersDimension();
/** Compute the partial derivatives of ODE with respect to state.
* @param t current value of the independent <I>time</I> variable
* @param y array containing the current value of the state vector
* @param yDot array containing the current value of the time derivative of the state vector
* @param dFdY placeholder array where to put the jacobian of the ODE with respect to the state vector
* @param dFdP placeholder array where to put the jacobian of the ODE with respect to the parameters
* @throws MathUserException this exception is propagated to the caller if the
* underlying user function triggers one
*/
void computeJacobians(double t, double[] y, double[] yDot, double[][] dFdY, double[][] dFdP)
throws MathUserException;
}

46
src/main/java/org/apache/commons/math/ode/jacobians/ParameterizedODE.java
View File

@ -1,46 +0,0 @@
/*
* 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.FirstOrderDifferentialEquations;
/** This interface represents {@link FirstOrderDifferentialEquations
* first order differential equations} with parameters.
*
* @see FirstOrderIntegratorWithJacobians
*
* @version $Revision$ $Date$
* @since 2.1
*/
public interface ParameterizedODE extends FirstOrderDifferentialEquations {
/** Get the number of parameters.
* @return number of parameters
*/
int getParametersDimension();
/** Set a parameter.
* @param i index of the parameters (must be between 0
* and {@link #getParametersDimension() getParametersDimension() - 1})
* @param value value for the parameter
*/
void setParameter(int i, double value);
}

95
src/main/java/org/apache/commons/math/ode/jacobians/StepHandlerWithJacobians.java
View File

@ -1,95 +0,0 @@
/*
* 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.exception.MathUserException;
/**
* This interface represents a handler that should be called after
* each successful step.
*
* <p>The ODE integrators compute the evolution of the state vector at
* some grid points that depend on their own internal algorithm. Once
* they have found a new grid point (possibly after having computed
* several evaluation of the derivative at intermediate points), they
* provide it to objects implementing this interface. These objects
* typically either ignore the intermediate steps and wait for the
* last one, store the points in an ephemeris, or forward them to
* specialized processing or output methods.</p>
*
* <p>Note that is is possible to register a {@link
* org.apache.commons.math.ode.sampling.StepHandler classical step handler}
* in the low level integrator used to build a {@link FirstOrderIntegratorWithJacobians}
* rather than implementing this class. The step handlers registered at low level
* will see the big compound state whether the step handlers defined by this interface
* see the original state, and its jacobians in separate arrays.</p>
*
* <p>The compound state is guaranteed to contain the original state in the first
* elements, followed by the jacobian with respect to initial state (in row order),
* followed by the jacobian with respect to parameters (in row order). If for example
* the original state dimension is 6 and there are 3 parameters, the compound state will
* be a 60 elements array. The first 6 elements will be the original state, the next 36
* elements will be the jacobian with respect to initial state, and the remaining 18 elements
* will be the jacobian with respect to parameters.</p>
*
* <p>Dealing with low level step handlers is cumbersome if one really needs the jacobians
* in these methods, but it also prevents many data being copied back and forth between
* state and jacobians on one side and compound state on the other side. So for performance
* reasons, it is recommended to use this interface <em>only</em> if jacobians are really
* needed and to use lower level handlers if only state is needed.</p>
*
* @see FirstOrderIntegratorWithJacobians
* @see StepInterpolatorWithJacobians
* @version $Revision$ $Date$
* @since 2.1
*/
public interface StepHandlerWithJacobians {
/** Determines whether this handler needs dense output.
* <p>This method allows the integrator to avoid performing extra
* computation if the handler does not need dense output.</p>
* @return true if the handler needs dense output
*/
boolean requiresDenseOutput();
/** Reset the step handler.
* Initialize the internal data as required before the first step is
* handled.
*/
void reset();
/**
* Handle the last accepted step
* @param interpolator interpolator for the last accepted step. For
* efficiency purposes, the various integrators reuse the same
* object on each call, so if the instance wants to keep it across
* all calls (for example to provide at the end of the integration a
* continuous model valid throughout the integration range, as the
* {@link org.apache.commons.math.ode.ContinuousOutputModel
* ContinuousOutputModel} class does), it should build a local copy
* using the clone method of the interpolator and store this copy.
* Keeping only a reference to the interpolator and reusing it will
* result in unpredictable behavior (potentially crashing the application).
* @param isLast true if the step is the last one
* @throws MathUserException this exception is propagated to the
* caller if the underlying user function triggers one
*/
void handleStep(StepInterpolatorWithJacobians interpolator, boolean isLast) throws MathUserException;
}

186
src/main/java/org/apache/commons/math/ode/jacobians/StepInterpolatorWithJacobians.java
View File

@ -1,186 +0,0 @@
/*
* 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 java.io.Externalizable;
import org.apache.commons.math.exception.MathUserException;
/** This interface represents an interpolator over the last step
* during an ODE integration.
*
* <p>The various ODE integrators provide objects implementing this
* interface to the step handlers. These objects are often custom
* objects tightly bound to the integrator internal algorithms. The
* handlers can use these objects to retrieve the state vector at
* intermediate times between the previous and the current grid points
* (this feature is often called dense output).</p>
* <p>One important thing to note is that the step handlers may be so
* tightly bound to the integrators that they often share some internal
* state arrays. This imply that one should <em>never</em> use a direct
* reference to a step interpolator outside of the step handler, either
* for future use or for use in another thread. If such a need arise, the
* step interpolator <em>must</em> be copied using the dedicated
* {@link #copy()} method.
* </p>
*
* @see FirstOrderIntegratorWithJacobians
* @see StepHandlerWithJacobians
* @version $Revision$ $Date$
* @since 2.1
*/
public interface StepInterpolatorWithJacobians extends Externalizable {
/**
* Get the previous grid point time.
* @return previous grid point time
*/
double getPreviousTime();
/**
* Get the current grid point time.
* @return current grid point time
*/
double getCurrentTime();
/**
* Get the time of the interpolated point.
* If {@link #setInterpolatedTime} has not been called, it returns
* the current grid point time.
* @return interpolation point time
*/
double getInterpolatedTime();
/**
* Set the time of the interpolated point.
* <p>Setting the time outside of the current step is now allowed, but
* should be used with care since the accuracy of the interpolator will
* probably be very poor far from this step. This allowance has been
* added to simplify implementation of search algorithms near the
* step endpoints.</p>
* <p>Setting the time changes the instance internal state. If a
* specific state must be preserved, a copy of the instance must be
* created using {@link #copy()}.</p>
* @param time time of the interpolated point
*/
void setInterpolatedTime(double time);
/**
* Get the state vector of the interpolated point.
* <p>The returned vector is a reference to a reused array, so
* it should not be modified and it should be copied if it needs
* to be preserved across several calls.</p>
* @return state vector at time {@link #getInterpolatedTime}
* @see #getInterpolatedYDot()
* @throws MathUserException if this call induces an automatic
* step finalization that throws one
*/
double[] getInterpolatedY() throws MathUserException;
/**
* Get the partial derivatives of the state vector with respect to
* the initial state of the interpolated point.
* <p>The returned vector is a reference to a reused array, so
* it should not be modified and it should be copied if it needs
* to be preserved across several calls.</p>
* @return partial derivatives of the state vector with respect to
* the initial state at time {@link #getInterpolatedTime}
* @see #getInterpolatedY()
* @throws MathUserException if this call induces an automatic
* step finalization that throws one
*/
double[][] getInterpolatedDyDy0() throws MathUserException;
/**
* Get the partial derivatives of the state vector with respect to
* the ODE parameters of the interpolated point.
* <p>The returned vector is a reference to a reused array, so
* it should not be modified and it should be copied if it needs
* to be preserved across several calls.</p>
* @return partial derivatives of the state vector with respect to
* the ODE parameters at time {@link #getInterpolatedTime}
* @see #getInterpolatedY()
* @throws MathUserException if this call induces an automatic
* step finalization that throws one
*/
double[][] getInterpolatedDyDp() throws MathUserException;
/**
* Get the time derivatives of the state vector of the interpolated point.
* <p>The returned vector is a reference to a reused array, so
* it should not be modified and it should be copied if it needs
* to be preserved across several calls.</p>
* @return derivatives of the state vector at time {@link #getInterpolatedTime}
* @see #getInterpolatedY()
* @throws MathUserException if this call induces an automatic
* step finalization that throws one
*/
double[] getInterpolatedYDot() throws MathUserException;
/**
* Get the time derivatives of the jacobian of the state vector
* with respect to the initial state of the interpolated point.
* <p>The returned vector is a reference to a reused array, so
* it should not be modified and it should be copied if it needs
* to be preserved across several calls.</p>
* @return time derivatives of the jacobian of the state vector
* with respect to the initial state at time {@link #getInterpolatedTime}
* @see #getInterpolatedY()
* @throws MathUserException if this call induces an automatic
* step finalization that throws one
*/
double[][] getInterpolatedDyDy0Dot() throws MathUserException;
/**
* Get the time derivatives of the jacobian of the state vector
* with respect to the ODE parameters of the interpolated point.
* <p>The returned vector is a reference to a reused array, so
* it should not be modified and it should be copied if it needs
* to be preserved across several calls.</p>
* @return time derivatives of the jacobian of the state vector
* with respect to the ODE parameters at time {@link #getInterpolatedTime}
* @see #getInterpolatedY()
* @throws MathUserException if this call induces an automatic
* step finalization that throws one
*/
double[][] getInterpolatedDyDpDot() throws MathUserException;
/** Check if the natural integration direction is forward.
* <p>This method provides the integration direction as specified by
* the integrator itself, it avoid some nasty problems in
* degenerated cases like null steps due to cancellation at step
* initialization, step control or discrete events
* triggering.</p>
* @return true if the integration variable (time) increases during
* integration
*/
boolean isForward();
/** Copy the instance.
* <p>The copied instance is guaranteed to be independent from the
* original one. Both can be used with different settings for
* interpolated time without any side effect.</p>
* @return a deep copy of the instance, which can be used independently.
* @throws MathUserException if this call induces an automatic
* step finalization that throws one
* @see #setInterpolatedTime(double)
*/
StepInterpolatorWithJacobians copy() throws MathUserException;
}

237
src/main/java/org/apache/commons/math/ode/jacobians/package.html
View File

@ -1,237 +0,0 @@
<html>
<!--
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.
-->
<!-- $Revision$ -->
<body>
<p>
This package provides classes to solve Ordinary Differential Equations problems
and also compute derivatives of the solution.
</p>
<p>
This package solves Initial Value Problems of the form
<code>y'=f(t,y,p)</code> with <code>t<sub>0</sub></code>,
<code>y(t<sub>0</sub>)=y<sub>0</sub></code> and
<code>dy(t<sub>0</sub>)/dp</sub></code> known. The provided
integrator computes estimates of <code>y(t)</code>,
<code>dy(t)/dy<sub>0</sub></code> and <code>dy(t)/dp</code>
from <code>t=t<sub>0</sub></code> to <code>t=t<sub>1</sub></code>,
where <code>y</code> is the state and <code>p</code> is a parameters
array.
</p>
<p>
The classes in this package mimic the behavior of classes and interfaces from the
<a href="../package-summary.html">org.apache.commons.math.ode</a>,
<a href="../events/package-summary.html">org.apache.commons.math.ode.events</a>
and <a href="../sampling/package-summary.html">org.apache.commons.math.ode.sampling</a>
packages, adding the jacobians <code>dy(t)/dy<sub>0</sub></code> and
<code>dy(t)/dp</code> to the methods signatures.
</p>
<p>
The classes and interfaces in this package mimic the behavior of the classes and
interfaces of the top level ode package, only adding parameters arrays for the jacobians.
The behavior of these classes is to create a compound state vector z containing both
the state y(t) and its derivatives dy(t)/dy<sub>0</sub> and dy(t<sub>0</sub>)/dp and
to set up an extended problem by adding the equations for the jacobians automatically.
These extended state and problems are then provided to a classical underlying integrator
chosen by user.
</p>
<p>
This behavior imply there will be a top level integrator knowing about state and jacobians
and a low level integrator knowing only about compound state (which may be big). If the user
wants to deal with the top level only, he will use the specialized step handler and event
handler classes registered at top level. He can also register classical step handlers and
event handlers, but in this case will see the big compound state. This state is guaranteed
to contain the original state in the first elements, followed by the jacobian with respect
to initial state (in row order), followed by the jacobian with respect to parameters (in
row order). If for example the original state dimension is 6 and there are 3 parameters,
the compound state will be a 60 elements array. The first 6 elements will be the original
state, the next 36 elements will be the jacobian with respect to initial state, and the
remaining 18 will be the jacobian with respect to parameters. Dealing with low level
step handlers and event handlers is cumbersome if one really needs the jacobians in these
methods, but it also prevents many data being copied back and forth between state and
jacobians on one side and compound state on the other side.
</p>
<p>
Here is a simple example of usage. We consider a two-dimensional problem where the
state vector y is the solution of the ordinary differential equations
<ul>
<li>y'<sub>0</sub>(t) = &omega; &times; (c<sub>1</sub> - y<sub>1</sub>(t))</li>
<li>y'<sub>1</sub>(t) = &omega; &times; (y<sub>0</sub>(t) - c<sub>0</sub>)</li>
</ul>
with some initial state y(t<sub>0</sub>) = (y<sub>0</sub>(t<sub>0</sub>),
y<sub>1</sub>(t<sub>0</sub>)).
</p>
<p>
The point trajectory depends on the initial state y(t<sub>0</sub>) and on the ODE
parameter &omega;. We want to compute both the final point position y(t<sub>end</sub>)
and the sensitivity of this point with respect to the initial state:
dy(t<sub>end</sub>)/dy(t<sub>0</sub>) which is a 2&times;2 matrix and its sensitivity
with respect to the parameter: dy(t<sub>end</sub>)/d&omega; which is a 2&times;1 matrix.
</p>
<p>
We consider first the simplest implementation: we define only the ODE and let
the classes compute the necessary jacobians by itself:
<code><pre>
public class BasicCircleODE implements ParameterizedODE {
private double[] c;
private double omega;
public BasicCircleODE(double[] c, double omega) {
this.c = c;
this.omega = omega;
}
public int getDimension() {
return 2;
}
public void computeDerivatives(double t, double[] y, double[] yDot) {
yDot[0] = omega * (c[1] - y[1]);
yDot[1] = omega * (y[0] - c[0]);
}
public int getParametersDimension() {
// we are only interested in the omega parameter
return 1;
}
public void setParameter(int i, double value) {
omega = value;
}
}
</pre></code>
</p>
<p>
We compute the results we want as follows:
<code><pre>
// low level integrator
FirstOrderIntegrator lowIntegrator =
new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10);
// set up ODE
double cx = 1.0;
double cy = 1.0;
double omega = 0.1;
ParameterizedODE ode = new BasicCircleODE(new double[] { cx, cy }, omega);
// set up high level integrator, using finite differences step hY and hP to compute jacobians
double[] hY = new double[] { 1.0e-5, 1.0e-5 };
double[] hP = new double[] { 1.0e-5 };
FirstOrderIntegratorWithJacobians integrator =
new FirstOrderIntegratorWithJacobians(lowIntegrator, ode, hY, hP);
// set up initial state and derivatives
double t0 = 0.0;
double[] y0 = new double[] { 0.0, cy };
double[][] dy0dp = new double[2][1] = { { 0.0 }, { 0.0 } }; // y0 does not depend on omega
// solve problem
double t = Math.PI / (2 * omega);
double[] y = new double[2];
double[][] dydy0 = new double[2][2];
double[][] dydp = new double[2][1];
integrator.integrate(t0, y0, dy0dp, t, y, dydy0, dydp);
</pre></code>
</p>
<p>
If in addition to getting the end state and its derivatives, we want to print the state
throughout integration process, we have to register a step handler. Inserting the following
before the call to integrate does the trick:
<code><pre>
StpeHandlerWithJacobians stepHandler = new StpeHandlerWithJacobians() {
public void reset() {}
public boolean requiresDenseOutput() { return false; }
public void handleStep(StepInterpolatorWithJacobians interpolator, boolean isLast)
throws DerivativeException {
double t = interpolator.getCurrentTime();
double[] y = interpolator.getInterpolatedY();
System.out.println(t + " " + y[0] + " " + y[1]);
}
};
integrator.addStepHandler(stepHandler);
</pre></code>
</p>
<p>
The implementation above relies on finite differences with small step sizes to compute the
internal jacobians. Since the ODE is really simple here, a better way is to compute them
exactly. So instead of implementing ParameterizedODE, we implement the ODEWithJacobians
interface as follows (i.e. we replace the setParameter method by a computeJacobians method):
<code><pre>
public class EnhancedCircleODE implements ODEWithJacobians {
private double[] c;
private double omega;
public EnhancedCircleODE(double[] c, double omega) {
this.c = c;
this.omega = omega;
}
public int getDimension() {
return 2;
}
public void computeDerivatives(double t, double[] y, double[] yDot) {
yDot[0] = omega * (c[1] - y[1]);
yDot[1] = omega * (y[0] - c[0]);
}
public int getParametersDimension() {
// we are only interested in the omega parameter
return 1;
}
public void computeJacobians(double t, double[] y, double[] yDot, double[][] dFdY, double[][] dFdP) {
dFdY[0][0] = 0;
dFdY[0][1] = -omega;
dFdY[1][0] = omega;
dFdY[1][1] = 0;
dFdP[0][0] = 0;
dFdP[0][1] = omega;
dFdP[0][2] = c[1] - y[1];
dFdP[1][0] = -omega;
dFdP[1][1] = 0;
dFdP[1][2] = y[0] - c[0];
}
}
</pre></code>
With this implementation, the hY and hP arrays are not needed anymore:
<code><pre>
ODEWithJacobians ode = new EnhancedCircleODE(new double[] { cx, cy }, omega);
FirstOrderIntegratorWithJacobians integrator =
new FirstOrderIntegratorWithJacobians(lowIntegrator, ode);
</pre></code>
</p>
</body>
</html>

8
src/site/xdoc/changes.xml
View File

@ -52,6 +52,9 @@ The <action> type attribute can be add,update,fix,remove.
If the output is not quite correct, check for invisible trailing spaces!
-->
<release version="3.0" date="TBD" description="TBD">
<action dev="luc" type="fix" issue="MATH-380">
Removed the ode.jacobians package.
</action>
<action dev="erans" type="update" issue="MATH-440">
Removed classes "FunctionEvaluationException", "MatrixVisitorException"
and "DerivativeException" (superseded by the new "MathUserException").
@ -105,6 +108,11 @@ The <action> type attribute can be add,update,fix,remove.
</action>
</release>
<release version="2.2" date="TBD" description="TBD">
<action dev="luc" type="fix" issue="MATH-380">
Deprecated the whole ode.jacobians package. It is clumsy and difficult to use. It will
be replaced by a completely rewritten implementation in 3.0, which will be more tightly
bound to the top level ode package
</action>
<action dev="luc" type="fix" issue="MATH-426" due-to="Erik van Ingen">
Added a normalization feature to transform samples so they have zero mean and unit standard deviation
</action>

436
src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java
View File

@ -1,436 +0,0 @@
/*
* 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.exception.MathUserException;
import org.apache.commons.math.ode.FirstOrderIntegrator;
import org.apache.commons.math.ode.IntegratorException;
import org.apache.commons.math.ode.nonstiff.DormandPrince54Integrator;
import org.apache.commons.math.stat.descriptive.SummaryStatistics;
import org.apache.commons.math.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
public class FirstOrderIntegratorWithJacobiansTest {
@Test
public void testLowAccuracyExternalDifferentiation()
throws IntegratorException, MathUserException {
// this test does not really test FirstOrderIntegratorWithJacobians,
// it only shows that WITHOUT this class, attempting to recover
// the jacobians from external differentiation on simple integration
// results with low accuracy gives very poor results. In fact,
// the curves dy/dp = g(b) when b varies from 2.88 to 3.08 are
// essentially noise.
// This test is taken from Hairer, Norsett and Wanner book
// Solving Ordinary Differential Equations I (Nonstiff problems),
// the curves dy/dp = g(b) are in figure 6.5
FirstOrderIntegrator integ =
new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 });
double hP = 1.0e-12;
SummaryStatistics residualsP0 = new SummaryStatistics();
SummaryStatistics residualsP1 = new SummaryStatistics();
for (double b = 2.88; b < 3.08; b += 0.001) {
Brusselator brusselator = new Brusselator(b);
double[] y = { 1.3, b };
integ.integrate(brusselator, 0, y, 20.0, y);
double[] yP = { 1.3, b + hP };
brusselator.setParameter(0, b + hP);
integ.integrate(brusselator, 0, yP, 20.0, yP);
residualsP0.addValue((yP[0] - y[0]) / hP - brusselator.dYdP0());
residualsP1.addValue((yP[1] - y[1]) / hP - brusselator.dYdP1());
}
Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) > 600);
Assert.assertTrue(residualsP0.getStandardDeviation() > 30);
Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) > 800);
Assert.assertTrue(residualsP1.getStandardDeviation() > 50);
}
@Test
public void testHighAccuracyExternalDifferentiation()
throws IntegratorException, MathUserException {
FirstOrderIntegrator integ =
new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 1.0e-10, 1.0e-10 });
double hP = 1.0e-12;
SummaryStatistics residualsP0 = new SummaryStatistics();
SummaryStatistics residualsP1 = new SummaryStatistics();
for (double b = 2.88; b < 3.08; b += 0.001) {
Brusselator brusselator = new Brusselator(b);
double[] y = { 1.3, b };
integ.integrate(brusselator, 0, y, 20.0, y);
double[] yP = { 1.3, b + hP };
brusselator.setParameter(0, b + hP);
integ.integrate(brusselator, 0, yP, 20.0, yP);
residualsP0.addValue((yP[0] - y[0]) / hP - brusselator.dYdP0());
residualsP1.addValue((yP[1] - y[1]) / hP - brusselator.dYdP1());
}
Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) > 0.02);
Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.03);
Assert.assertTrue(residualsP0.getStandardDeviation() > 0.003);
Assert.assertTrue(residualsP0.getStandardDeviation() < 0.004);
Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) > 0.04);
Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05);
Assert.assertTrue(residualsP1.getStandardDeviation() > 0.007);
Assert.assertTrue(residualsP1.getStandardDeviation() < 0.008);
}
@Test
public void testInternalDifferentiation()
throws IntegratorException, MathUserException {
FirstOrderIntegrator integ =
new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 });
double hP = 1.0e-12;
SummaryStatistics residualsP0 = new SummaryStatistics();
SummaryStatistics residualsP1 = new SummaryStatistics();
for (double b = 2.88; b < 3.08; b += 0.001) {
Brusselator brusselator = new Brusselator(b);
brusselator.setParameter(0, b);
double[] z = { 1.3, b };
double[][] dZdZ0 = new double[2][2];
double[][] dZdP = new double[2][1];
double hY = 1.0e-12;
FirstOrderIntegratorWithJacobians extInt =
new FirstOrderIntegratorWithJacobians(integ, brusselator, new double[] { b },
new double[] { hY, hY }, new double[] { hP });
extInt.setMaxEvaluations(5000);
extInt.integrate(0, z, new double[][] { { 0.0 }, { 1.0 } }, 20.0, z, dZdZ0, dZdP);
Assert.assertEquals(5000, extInt.getMaxEvaluations());
Assert.assertTrue(extInt.getEvaluations() > 1500);
Assert.assertTrue(extInt.getEvaluations() < 2100);
Assert.assertEquals(4 * integ.getEvaluations(), extInt.getEvaluations());
residualsP0.addValue(dZdP[0][0] - brusselator.dYdP0());
residualsP1.addValue(dZdP[1][0] - brusselator.dYdP1());
}
Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.02);
Assert.assertTrue(residualsP0.getStandardDeviation() < 0.003);
Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05);
Assert.assertTrue(residualsP1.getStandardDeviation() < 0.01);
}
@Test
public void testAnalyticalDifferentiation()
throws IntegratorException, MathUserException {
FirstOrderIntegrator integ =
new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 });
SummaryStatistics residualsP0 = new SummaryStatistics();
SummaryStatistics residualsP1 = new SummaryStatistics();
for (double b = 2.88; b < 3.08; b += 0.001) {
Brusselator brusselator = new Brusselator(b);
brusselator.setParameter(0, b);
double[] z = { 1.3, b };
double[][] dZdZ0 = new double[2][2];
double[][] dZdP = new double[2][1];
FirstOrderIntegratorWithJacobians extInt =
new FirstOrderIntegratorWithJacobians(integ, brusselator);
extInt.setMaxEvaluations(5000);
extInt.integrate(0, z, new double[][] { { 0.0 }, { 1.0 } }, 20.0, z, dZdZ0, dZdP);
Assert.assertEquals(5000, extInt.getMaxEvaluations());
Assert.assertTrue(extInt.getEvaluations() > 350);
Assert.assertTrue(extInt.getEvaluations() < 510);
Assert.assertEquals(integ.getEvaluations(), extInt.getEvaluations());
residualsP0.addValue(dZdP[0][0] - brusselator.dYdP0());
residualsP1.addValue(dZdP[1][0] - brusselator.dYdP1());
}
Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.014);
Assert.assertTrue(residualsP0.getStandardDeviation() < 0.003);
Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05);
Assert.assertTrue(residualsP1.getStandardDeviation() < 0.01);
}
@Test
public void testFinalResult() throws IntegratorException, MathUserException {
FirstOrderIntegrator integ =
new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 1.0e-10, 1.0e-10 });
double[] y = new double[] { 0.0, 1.0 };
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 * FastMath.PI;
FirstOrderIntegratorWithJacobians extInt =
new FirstOrderIntegratorWithJacobians(integ, circle);
extInt.integrate(0, y, circle.exactDyDp(0), t, y, dydy0, dydp);
for (int i = 0; i < y.length; ++i) {
Assert.assertEquals(circle.exactY(t)[i], y[i], 1.0e-9);
}
for (int i = 0; i < dydy0.length; ++i) {
for (int j = 0; j < dydy0[i].length; ++j) {
Assert.assertEquals(circle.exactDyDy0(t)[i][j], dydy0[i][j], 1.0e-9);
}
}
for (int i = 0; i < dydp.length; ++i) {
for (int j = 0; j < dydp[i].length; ++j) {
Assert.assertEquals(circle.exactDyDp(t)[i][j], dydp[i][j], 1.0e-7);
}
}
}
@Test
public void testStepHandlerResult() throws IntegratorException, MathUserException {
FirstOrderIntegrator integ =
new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 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 * FastMath.PI;
final FirstOrderIntegratorWithJacobians extInt =
new FirstOrderIntegratorWithJacobians(integ, circle);
extInt.addStepHandler(new StepHandlerWithJacobians() {
public void reset() {
}
public boolean requiresDenseOutput() {
return false;
}
public void handleStep(StepInterpolatorWithJacobians interpolator, boolean isLast)
throws MathUserException {
double t = interpolator.getCurrentTime();
double[] y = interpolator.getInterpolatedY();
double[][] dydy0 = interpolator.getInterpolatedDyDy0();
double[][] dydp = interpolator.getInterpolatedDyDp();
Assert.assertEquals(interpolator.getPreviousTime(), extInt.getCurrentStepStart(), 1.0e-10);
Assert.assertTrue(extInt.getCurrentSignedStepsize() < 0.5);
for (int i = 0; i < y.length; ++i) {
Assert.assertEquals(circle.exactY(t)[i], y[i], 1.0e-9);
}
for (int i = 0; i < dydy0.length; ++i) {
for (int j = 0; j < dydy0[i].length; ++j) {
Assert.assertEquals(circle.exactDyDy0(t)[i][j], dydy0[i][j], 1.0e-9);
}
}
for (int i = 0; i < dydp.length; ++i) {
for (int j = 0; j < dydp[i].length; ++j) {
Assert.assertEquals(circle.exactDyDp(t)[i][j], dydp[i][j], 3.0e-8);
}
}
double[] yDot = interpolator.getInterpolatedYDot();
double[][] dydy0Dot = interpolator.getInterpolatedDyDy0Dot();
double[][] dydpDot = interpolator.getInterpolatedDyDpDot();
for (int i = 0; i < yDot.length; ++i) {
Assert.assertEquals(circle.exactYDot(t)[i], yDot[i], 1.0e-10);
}
for (int i = 0; i < dydy0Dot.length; ++i) {
for (int j = 0; j < dydy0Dot[i].length; ++j) {
Assert.assertEquals(circle.exactDyDy0Dot(t)[i][j], dydy0Dot[i][j], 1.0e-10);
}
}
for (int i = 0; i < dydpDot.length; ++i) {
for (int j = 0; j < dydpDot[i].length; ++j) {
Assert.assertEquals(circle.exactDyDpDot(t)[i][j], dydpDot[i][j], 3.0e-9);
}
}
}
});
extInt.integrate(0, y, circle.exactDyDp(0), t, y, dydy0, dydp);
}
@Test
public void testEventHandler() throws IntegratorException, MathUserException {
FirstOrderIntegrator integ =
new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 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 * FastMath.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 * FastMath.PI);
}
private static class Brusselator implements ParameterizedODE, ODEWithJacobians {
private double b;
public Brusselator(double b) {
this.b = b;
}
public int getDimension() {
return 2;
}
public void setParameter(int i, double p) {
b = p;
}
public int getParametersDimension() {
return 1;
}
public void computeDerivatives(double t, double[] y, double[] yDot) {
double prod = y[0] * y[0] * y[1];
yDot[0] = 1 + prod - (b + 1) * y[0];
yDot[1] = b * y[0] - prod;
}
public void computeJacobians(double t, double[] y, double[] yDot, double[][] dFdY, double[][] dFdP) {
double p = 2 * y[0] * y[1];
double y02 = y[0] * y[0];
dFdY[0][0] = p - (1 + b);
dFdY[0][1] = y02;
dFdY[1][0] = b - p;
dFdY[1][1] = -y02;
dFdP[0][0] = -y[0];
dFdP[1][0] = y[0];
}
public double dYdP0() {
return -1088.232716447743 + (1050.775747149553 + (-339.012934631828 + 36.52917025056327 * b) * b) * b;
}
public double dYdP1() {
return 1502.824469929139 + (-1438.6974831849952 + (460.959476642384 - 49.43847385647082 * b) * b) * b;
}
}
/** ODE representing a point moving on a circle with provided center and angular rate. */
private static class Circle implements ODEWithJacobians {
private final double[] y0;
private double cx;
private double cy;
private double omega;
public Circle(double[] y0, double cx, double cy, double omega) {
this.y0 = y0.clone();
this.cx = cx;
this.cy = cy;
this.omega = omega;
}
public int getDimension() {
return 2;
}
public int getParametersDimension() {
return 3;
}
public void computeDerivatives(double t, double[] y, double[] yDot) {
yDot[0] = omega * (cy - y[1]);
yDot[1] = omega * (y[0] - cx);
}
public void computeJacobians(double t, double[] y, double[] yDot, double[][] dFdY, double[][] dFdP) {
dFdY[0][0] = 0;
dFdY[0][1] = -omega;
dFdY[1][0] = omega;
dFdY[1][1] = 0;
dFdP[0][0] = 0;
dFdP[0][1] = omega;
dFdP[0][2] = cy - y[1];
dFdP[1][0] = -omega;
dFdP[1][1] = 0;
dFdP[1][2] = y[0] - cx;
}
public double[] exactY(double t) {
double cos = FastMath.cos(omega * t);
double sin = FastMath.sin(omega * t);
double dx0 = y0[0] - cx;
double dy0 = y0[1] - cy;
return new double[] {
cx + cos * dx0 - sin * dy0,
cy + sin * dx0 + cos * dy0
};
}
public double[][] exactDyDy0(double t) {
double cos = FastMath.cos(omega * t);
double sin = FastMath.sin(omega * t);
return new double[][] {
{ cos, -sin },
{ sin, cos }
};
}
public double[][] exactDyDp(double t) {
double cos = FastMath.cos(omega * t);
double sin = FastMath.sin(omega * t);
double dx0 = y0[0] - cx;
double dy0 = y0[1] - cy;
return new double[][] {
{ 1 - cos, sin, -t * (sin * dx0 + cos * dy0) },
{ -sin, 1 - cos, t * (cos * dx0 - sin * dy0) }
};
}
public double[] exactYDot(double t) {
double oCos = omega * FastMath.cos(omega * t);
double oSin = omega * FastMath.sin(omega * t);
double dx0 = y0[0] - cx;
double dy0 = y0[1] - cy;
return new double[] {
-oSin * dx0 - oCos * dy0,
oCos * dx0 - oSin * dy0
};
}
public double[][] exactDyDy0Dot(double t) {
double oCos = omega * FastMath.cos(omega * t);
double oSin = omega * FastMath.sin(omega * t);
return new double[][] {
{ -oSin, -oCos },
{ oCos, -oSin }
};
}
public double[][] exactDyDpDot(double t) {
double cos = FastMath.cos(omega * t);
double sin = FastMath.sin(omega * t);
double oCos = omega * cos;
double oSin = omega * sin;
double dx0 = y0[0] - cx;
double dy0 = y0[1] - cy;
return new double[][] {
{ oSin, oCos, -sin * dx0 - cos * dy0 - t * ( oCos * dx0 - oSin * dy0) },
{ -oCos, oSin, cos * dx0 - sin * dy0 + t * (-oSin * dx0 - oCos * dy0) }
};
}
}
}