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

added support for secondary state in ODE step interpolators 

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1176734 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2011-09-28 05:56:42 +00:00
parent 924769fb05
commit b5724e1cc0
32 changed files with 1228 additions and 1228 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java
View File

@ -315,9 +315,8 @@ public enum LocalizedFormats implements Localizable {
UNABLE_TO_SOLVE_SINGULAR_PROBLEM("unable to solve: singular problem"),
UNBOUNDED_SOLUTION("unbounded solution"),
UNKNOWN_MODE("unknown mode {0}, known modes: {1} ({2}), {3} ({4}), {5} ({6}), {7} ({8}), {9} ({10}) and {11} ({12})"),
UNKNOWN_ADDITIONAL_EQUATION("unknown additional equation"),
UNKNOWN_PARAMETER("unknown parameter {0}"),
UNMATCHED_ODE_IN_EXTENDED_SET("ode does not match the main ode set in the extended set"),
UNMATCHED_ODE_IN_EXPANDED_SET("ode does not match the main ode set in the extended set"),
CANNOT_PARSE_AS_TYPE("string {0} unparseable (from position {1}) as an object of type {2}"), /* keep */
CANNOT_PARSE("string {0} unparseable (from position {1})"), /* keep */
UNPARSEABLE_3D_VECTOR("unparseable 3D vector: \"{0}\""),

86
src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java
View File

@ -47,7 +47,7 @@ import org.apache.commons.math.util.MathUtils;
* @version $Id$
* @since 2.0
*/
public abstract class AbstractIntegrator implements ExpandableFirstOrderIntegrator {
public abstract class AbstractIntegrator implements FirstOrderIntegrator {
/** Step handler. */
protected Collection<StepHandler> stepHandlers;
@ -77,7 +77,7 @@ public abstract class AbstractIntegrator implements ExpandableFirstOrderIntegrat
private Incrementor evaluations;
/** Differential equations to integrate. */
private transient ExpandableFirstOrderDifferentialEquations equations;
private transient ExpandableStatefulODE equations;
/** Build an instance.
* @param name name of the method
@ -185,21 +185,59 @@ public abstract class AbstractIntegrator implements ExpandableFirstOrderIntegrat
evaluations.resetCount();
}
/** Set the differential equations.
* @param equations differential equations to integrate
* @see #computeDerivatives(double, double[], double[])
/** Set the equations.
* @param equations equations to set
*/
protected void setEquations(final ExpandableFirstOrderDifferentialEquations equations) {
protected void setEquations(final ExpandableStatefulODE equations) {
this.equations = equations;
}
/** {@inheritDoc} */
public double integrate(FirstOrderDifferentialEquations equations,
double t0, double[] y0, double t, double[] y)
public double integrate(final FirstOrderDifferentialEquations equations,
final double t0, final double[] y0, final double t, final double[] y)
throws MathIllegalStateException, MathIllegalArgumentException {
return integrate(new ExpandableFirstOrderDifferentialEquations(equations), t0, y0, t, y);
if (y0.length != equations.getDimension()) {
throw new DimensionMismatchException(y0.length, equations.getDimension());
}
if (y.length != equations.getDimension()) {
throw new DimensionMismatchException(y.length, equations.getDimension());
}
// prepare expandable stateful equations
final ExpandableStatefulODE expandable = new ExpandableStatefulODE(equations);
expandable.setTime(t0);
expandable.setPrimaryState(y0);
// perform integration
integrate(expandable, t);
// extract results back from the stateful equations
System.arraycopy(expandable.getPrimaryState(), 0, y, 0, y.length);
return expandable.getTime();
}
/** Integrate a set of differential equations up to the given time.
* <p>This method solves an Initial Value Problem (IVP).</p>
* <p>The set of differential equations is composed of a main set, which
* can be extended by some sets of secondary equations. The set of
* equations must be already set up with initial time and partial states.
* At integration completion, the final time and partial states will be
* available in the same object.</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 equations complete set of differential equations to integrate
* @param t target time for the integration
* (can be set to a value smaller than <code>t0</code> for backward integration)
* @throws MathIllegalStateException if the integrator cannot perform integration
* @throws MathIllegalArgumentException if integration parameters are wrong (typically
* too small integration span)
*/
public abstract void integrate(ExpandableStatefulODE equations, double t)
throws MathIllegalStateException, MathIllegalArgumentException;
/** Compute the derivatives and check the number of evaluations.
* @param t current value of the independent <I>time</I> variable
* @param y array containing the current value of the state vector
@ -335,33 +373,19 @@ public abstract class AbstractIntegrator implements ExpandableFirstOrderIntegrat
}
/** Perform some sanity checks on the integration parameters.
* @param ode differential equations set
* @param t0 start time
* @param y0 state vector at t0
/** Check the integration span.
* @param t target time for the integration
* @param y placeholder where to put the state vector
* @exception DimensionMismatchException if some inconsistency is detected
* @exception NumberIsTooSmallException if integration span is too small
*/
protected void sanityChecks(final ExpandableFirstOrderDifferentialEquations ode,
final double t0, final double[] y0,
final double t, final double[] y)
throws DimensionMismatchException, NumberIsTooSmallException {
protected void sanityChecks(final ExpandableStatefulODE equations, final double t)
throws NumberIsTooSmallException {
if (ode.getMainSetDimension() != y0.length) {
throw new DimensionMismatchException(ode.getDimension(), y0.length);
}
if (ode.getMainSetDimension() != y.length) {
throw new DimensionMismatchException(ode.getDimension(), y.length);
}
if (FastMath.abs(t - t0) <= 1.0e-12 * FastMath.max(FastMath.abs(t0), FastMath.abs(t))) {
final double threshold = 1000 * FastMath.ulp(FastMath.max(FastMath.abs(equations.getTime()),
FastMath.abs(t)));
final double dt = FastMath.abs(equations.getTime() - t);
if (dt <= threshold) {
throw new NumberIsTooSmallException(LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL,
FastMath.abs(t - t0),
1.0e-12 * FastMath.max(FastMath.abs(t0), FastMath.abs(t)),
false);
dt, threshold, false);
}
}

88
src/main/java/org/apache/commons/math/ode/AdditionalStateAndEquations.java
View File

@ -1,88 +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;
/**
* This class is a container for additional state parameters and their associated
* evolution equation.
* <p>
* It is a container allowing the integrator to keep constant consistency between
* additional states and the corresponding equations. It allows to set additional
* state values, get current additional state value and derivatives by reference
* on the associated additional equations.
* </p>
*
* @see ExpandableFirstOrderDifferentialEquations
* @see AdditionalEquations
*
* @version $Id$
* @since 3.0
*/
class AdditionalStateAndEquations {
/** Additional equations set. */
private final AdditionalEquations addEquations;
/** Current additional state. */
private double[] addState;
/** Current additional state derivatives. */
private double[] addStateDot;
/** Create a new instance based on one set of additional equations.
* @param addEqu additional equations.
*/
public AdditionalStateAndEquations(final AdditionalEquations addEqu) {
this.addEquations = addEqu;
}
/** Get a reference to the current value of the additional state.
* <p>The array returned is a true reference to the state array, so it may be
* used to store data into it.</>
* @return a reference current value of the additional state.
*/
public double[] getAdditionalState() {
return addState;
}
/** Get a reference to the current value of the additional state derivatives.
* <p>The array returned is a true reference to the state array, so it may be
* used to store data into it.</>
* @return a reference current value of the additional state derivatives.
*/
public double[] getAdditionalStateDot() {
return addStateDot;
}
/** Get the instance of the current additional equations.
* @return current value of the additional equations.
*/
public AdditionalEquations getAdditionalEquations() {
return addEquations;
}
/** Set a value to additional state.
* @param state additional state value.
*/
public void setAdditionalState(final double[] state) {
this.addState = state.clone();
this.addStateDot = new double[state.length];
}
}

98
src/main/java/org/apache/commons/math/ode/EquationsMapper.java Normal file
View File

@ -0,0 +1,98 @@
/*
* 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;
import java.io.Serializable;
import org.apache.commons.math.exception.DimensionMismatchException;
/**
* Class mapping the part of a complete state or derivative that pertains
* to a specific differential equation.
* <p>
* Instances of this class are guaranteed to be immutable.
* </p>
* @see SecondaryEquations
* @version $Id$
* @since 3.0
*/
public class EquationsMapper implements Serializable {
/** Serializable UID. */
private static final long serialVersionUID = 20110925L;
/** Index of the first equation element in complete state arrays. */
final int firstIndex;
/** Dimension of the secondary state parameters. */
final int dimension;
/** simple constructor.
* @param firstIndex index of the first equation element in complete state arrays
* @param dimension dimension of the secondary state parameters
*/
public EquationsMapper(final int firstIndex, final int dimension) {
this.firstIndex = firstIndex;
this.dimension = dimension;
}
/** Get the index of the first equation element in complete state arrays.
* @return index of the first equation element in complete state arrays
*/
public int getFirstIndex() {
return firstIndex;
}
/** Get the dimension of the secondary state parameters.
* @return dimension of the secondary state parameters
*/
public int getDimension() {
return dimension;
}
/** Extract equation data from a complete state or derivative array.
* @param complete complete state or derivative array from which
* equation data should be retrieved
* @param equationData placeholder where to put equation data
* @throws DimensionMismatchException if the dimension of the equation data does not
* match the mapper dimension
*/
public void extractEquationData(double[] complete, double[] equationData)
throws DimensionMismatchException {
if (equationData.length != dimension) {
throw new DimensionMismatchException(equationData.length, dimension);
}
System.arraycopy(complete, firstIndex, equationData, 0, dimension);
}
/** Insert equation data into a complete state or derivative array.
* @param equationData equation data to be inserted into the complete array
* @param complete placeholder where to put equation data (only the
* part corresponding to the equation will be overwritten)
* @throws DimensionMismatchException if the dimension of the equation data does not
* match the mapper dimension
*/
public void insertEquationData(double[] equationData, double[] complete)
throws DimensionMismatchException {
if (equationData.length != dimension) {
throw new DimensionMismatchException(equationData.length, dimension);
}
System.arraycopy(equationData, 0, complete, firstIndex, dimension);
}
}

272
src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderDifferentialEquations.java
View File

@ -1,272 +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;
import java.util.ArrayList;
import java.util.List;
import org.apache.commons.math.exception.MathIllegalArgumentException;
import org.apache.commons.math.exception.util.LocalizedFormats;
/**
* This class represents a combined set of first order differential equations,
* with at least a main set of equations expandable by some sets of additional
* equations.
* <p>
* This class extends the {@link FirstOrderDifferentialEquations}. It allows to
* identify which part of a complete set of differential equations correspond to
* the main set and which part correspond to the expansion sets.
* </p>
* <p>
* One typical use case is the computation of the jacobian matrix for some ODE.
* The main set of equations corresponds to the raw ODE, and we add to this set
* another bunch of equations which represent the jacobian matrix of the main
* set. In that case, we want the integrator to use <em>only</em> the main set
* to estimate the errors and hence the step sizes. It should <em>not</em> use
* the additional equations in this computation.
* The {@link ExpandableFirstOrderIntegrator integrator} will be able to know
* where the main set ends and so where the expansion sets begin.
* </p>
* <p>
* We consider that the main set always corresponds to the first equations and
* the expansion sets to the last equations.
* </p>
*
* @see FirstOrderDifferentialEquations
* @see JacobianMatrices
*
* @version $Id$
* @since 3.0
*/
public class ExpandableFirstOrderDifferentialEquations implements FirstOrderDifferentialEquations {
/** Main set of differential equations. */
private final FirstOrderDifferentialEquations mainSet;
/** Additional sets of equations and associated states. */
private final List<AdditionalStateAndEquations> addedSets;
/** Create a new instance of ExpandableFirstOrderDifferentialEquations.
* @param fode the main set of differential equations to be integrated.
*/
public ExpandableFirstOrderDifferentialEquations(final FirstOrderDifferentialEquations fode) {
this.mainSet = fode;
this.addedSets = new ArrayList<AdditionalStateAndEquations>();
}
/** Return the dimension of the whole state vector.
* <p>
* The whole state vector results in the assembly of the main set of
* equations and, if there are some, the added sets of equations.
* </p>
* @return dimension of the whole state vector
*/
public int getDimension()
throws MathIllegalArgumentException {
int dimension = this.getMainSetDimension();
try {
for (AdditionalStateAndEquations stateAndEqu : addedSets) {
dimension += stateAndEqu.getAdditionalEquations().getDimension();
}
return dimension;
} catch (Exception e) {
// TODO we should not catch Exception, and we should identify the offending additional equation
throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_ADDITIONAL_EQUATION);
}
}
/** Return the dimension of the main set of equations.
* <p>
* The main set of equations represents the first part of an ODE state.
* The error estimations and adaptive step size computation should be
* done on this first part only, not on the final part of the state
* which represents expansion sets of equations considered as secondary.
* </p>
* @return dimension of the main set of equations, must be lesser than or
* equal to the {@link #getDimension() total dimension}
*/
public int getMainSetDimension() {
return mainSet.getDimension();
}
/** Return the cumulated dimension of all added sets of equations.
* @return dimension of all added sets of equations
* @throws IllegalArgumentException if some additional equation is unknown
*/
public int getAddedSetsDimension()
throws IllegalArgumentException {
int addDim = 0;
try {
for (AdditionalStateAndEquations stateAndEqu : addedSets) {
addDim += stateAndEqu.getAdditionalEquations().getDimension();
}
return addDim;
} catch (Exception e) {
// TODO we should not catch Exception, and we should identify the offending additional equation
throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_ADDITIONAL_EQUATION);
}
}
/** Return the dimension of one added set of equations.
* @param addEqu Additional equations used as a reference for selection
* @return dimension of the added set of equations
* @throws IllegalArgumentException if additional equation is unknown
*/
public int getAddedSetDimension(final AdditionalEquations addEqu) {
return selectStateAndEquations(addEqu).getAdditionalEquations().getDimension();
}
/** Get the current time derivative of the total state vector.
* @param t current value of the independent <I>time</I> variable
* @param y array containing the current value of the state vector
* @param yDot placeholder array where to put the time derivative of the state vector
*/
public void computeDerivatives(final double t, final double[] y, final double[] yDot) {
// Add contribution for the main set of equations
int index = getMainSetDimension();
double[] m = new double[index];
double[] mDot = new double[index];
// update current main state
System.arraycopy(y, 0, m, 0, index);
// compute derivatives
mainSet.computeDerivatives(t, m, mDot);
// update main state contribution in global array
System.arraycopy(mDot, 0, yDot, 0, index);
// Add contribution for additional equations
for (final AdditionalStateAndEquations stateAndEqu : addedSets) {
final double[] p = stateAndEqu.getAdditionalState();
final double[] pDot = stateAndEqu.getAdditionalStateDot();
// update current additional state
System.arraycopy(y, index, p, 0, p.length);
// compute additional derivatives
stateAndEqu.getAdditionalEquations().computeDerivatives(t, m, mDot, p, pDot);
// update each additional state contribution in global array
System.arraycopy(pDot, 0, yDot, index, p.length);
// incrementing index
index += p.length;
}
}
/** Add a set of user-specified equations to be integrated along with
* the main set of equations.
*
* @param addEqu additional equations
* @see #setInitialAdditionalState(double[], AdditionalEquations)
* @see #getCurrentAdditionalState(AdditionalEquations)
*/
public void addAdditionalEquations(final AdditionalEquations addEqu) {
addedSets.add(new AdditionalStateAndEquations(addEqu));
}
/** Get the instance of the main set of equations.
* @return current value of the main set of equations.
*/
public FirstOrderDifferentialEquations getMainSet() {
return mainSet;
}
/** Set initial additional state.
* @param addState additional state
* @param addEqu additional equations used as a reference for selection
* @throws IllegalArgumentException if additional equation is unknown
*/
public void setInitialAdditionalState(final double[] addState, final AdditionalEquations addEqu) {
selectStateAndEquations(addEqu).setAdditionalState(addState);
}
/** Set current additional state.
* <p>
* The total current state computed by the integrator
* is dispatched here to the various additional states.
* </p>
* @param currentState total current state
* @throws IllegalArgumentException if additional equation is unknown
*/
public void setCurrentAdditionalState(final double[] currentState)
throws IllegalArgumentException {
int index = getMainSetDimension();
try {
for (AdditionalStateAndEquations stateAndEqu : addedSets) {
final int addDim = stateAndEqu.getAdditionalEquations().getDimension();
final double[] addState = new double[addDim];
System.arraycopy(currentState, index, addState, 0, addDim);
stateAndEqu.setAdditionalState(addState);
index += addDim;
}
} catch (Exception e) {
// TODO we should not catch Exception, and we should identify the offending additional equation
throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_ADDITIONAL_EQUATION);
}
}
/** Get current additional state.
* @param addEqu additional equations used as a reference for selection
* @return current additional state
* @throws IllegalArgumentException if additional equation is unknown
*/
public double[] getCurrentAdditionalState(final AdditionalEquations addEqu) {
return selectStateAndEquations(addEqu).getAdditionalState();
}
/** Get all current additional states accumulated.
* @return current additional states
* @throws IllegalArgumentException if additional equation is unknown
*/
public double[] getCurrentAdditionalStates()
throws IllegalArgumentException {
int index = 0;
final double[] cumulState = new double[getAddedSetsDimension()];
try {
for (AdditionalStateAndEquations stateAndEqu : addedSets) {
final int addDim = stateAndEqu.getAdditionalEquations().getDimension();
final double[] addState = stateAndEqu.getAdditionalState();
System.arraycopy(addState, 0, cumulState, index, addDim);
index += addDim;
}
return cumulState;
} catch (Exception e) {
// TODO we should not catch Exception, and we should identify the offending additional equation
throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_ADDITIONAL_EQUATION);
}
}
/** Select additional state and equations pair in the list.
* @param addEqu Additional equations used as a reference for selection
* @return additional state and equations pair
* @throws IllegalArgumentException if additional equation is unknown
*/
private AdditionalStateAndEquations selectStateAndEquations(final AdditionalEquations addEqu)
throws IllegalArgumentException {
for (AdditionalStateAndEquations stateAndEqu : addedSets) {
if (stateAndEqu.getAdditionalEquations() == addEqu) {
return stateAndEqu;
}
}
// TODO we should not catch Exception, and we should identify the offending additional equation
throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_ADDITIONAL_EQUATION);
}
}

65
src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderIntegrator.java
View File

@ -1,65 +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;
import org.apache.commons.math.exception.MathIllegalArgumentException;
import org.apache.commons.math.exception.MathIllegalStateException;
/**
* This interface represents a first order integrator for expandable
* differential equations.
* <p>
* The classes devoted to solve expandable first order differential equations
* should implement this interface. The problems which can be handled should
* implement the {@link ExpandableFirstOrderDifferentialEquations} interface.
* </p>
*
* @see ExpandableFirstOrderDifferentialEquations
*
* @version $Id$
* @since 3.0
*/
public interface ExpandableFirstOrderIntegrator extends FirstOrderIntegrator {
/** Integrate a set of differential equations up to the given time.
* <p>This method solves an Initial Value Problem (IVP).</p>
* <p>The set of differential equations is composed of a main set, which
* can be extended by some sets of additional equations.</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 equations complete set of differential equations to integrate
* @param t0 initial time
* @param y0 initial value of the main state vector 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 main state vector at each successful
* step (and hence at the end of integration), can be the same object as y0
* @return stop time, will be the same as target time if integration reached its
* target, but may be different if some {@link
* org.apache.commons.math.ode.events.EventHandler} stops it at some point.
* @throws MathIllegalStateException if the integrator cannot perform integration
* @throws MathIllegalArgumentException if integration parameters are wrong (typically
* too small integration span)
*/
double integrate(ExpandableFirstOrderDifferentialEquations equations,
double t0, double[] y0, double t, double[] y)
throws MathIllegalStateException, MathIllegalArgumentException;
}

327
src/main/java/org/apache/commons/math/ode/ExpandableStatefulODE.java Normal file
View File

@ -0,0 +1,327 @@
/*
* 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;
import java.util.ArrayList;
import java.util.List;
import org.apache.commons.math.exception.DimensionMismatchException;
/**
* This class represents a combined set of first order differential equations,
* with at least a primary set of equations expandable by some sets of secondary
* equations.
* <p>
* One typical use case is the computation of the Jacobian matrix for some ODE.
* In this case, the primary set of equations corresponds to the raw ODE, and we
* add to this set another bunch of secondary equations which represent the Jacobian
* matrix of the primary set.
* </p>
* <p>
* We want the integrator to use <em>only</em> the primary set to estimate the
* errors and hence the step sizes. It should <em>not</em> use the secondary
* equations in this computation. The {@link AbstractIntegrator integrator} will
* be able to know where the primary set ends and so where the secondary sets begin.
* </p>
*
* @see FirstOrderDifferentialEquations
* @see JacobianMatrices
*
* @version $Id$
* @since 3.0
*/
public class ExpandableStatefulODE {
/** Primary differential equation. */
private final FirstOrderDifferentialEquations primary;
/** Mapper for primary equation. */
private final EquationsMapper primaryMapper;
/** Time. */
private double time;
/** State. */
private final double[] primaryState;
/** State derivative. */
private final double[] primaryStateDot;
/** Components of the expandable ODE. */
private List<SecondaryComponent> components;
/** Build an expandable set from its primary ODE set.
* @param ode the primary set of differential equations to be integrated.
*/
public ExpandableStatefulODE(final FirstOrderDifferentialEquations primary) {
final int n = primary.getDimension();
this.primary = primary;
this.primaryMapper = new EquationsMapper(0, n);
this.time = Double.NaN;
this.primaryState = new double[n];
this.primaryStateDot = new double[n];
this.components = new ArrayList<ExpandableStatefulODE.SecondaryComponent>();
}
/** Get the primary set of differential equations.
* @return primary set of differential equations
*/
public FirstOrderDifferentialEquations getPrimary() {
return primary;
}
/** Return the dimension of the complete set of equations.
* <p>
* The complete set of equations correspond to the primary set plus all secondary sets.
* </p>
* @return dimension of the complete set of equations
*/
public int getTotalDimension() {
if (components.isEmpty()) {
// there are no secondary equations, the complete set is limited to the primary set
return primaryMapper.getDimension();
} else {
// there are secondary equations, the complete set ends after the last set
final EquationsMapper lastMapper = components.get(components.size() - 1).mapper;
return lastMapper.getFirstIndex() + lastMapper.getDimension();
}
}
/** Get the current time derivative of the complete state vector.
* @param t current value of the independent <I>time</I> variable
* @param y array containing the current value of the complete state vector
* @param yDot placeholder array where to put the time derivative of the complete state vector
*/
public void computeDerivatives(final double t, final double[] y, final double[] yDot) {
// compute derivatives of the primary equations
primaryMapper.extractEquationData(y, primaryState);
primary.computeDerivatives(t, primaryState, primaryStateDot);
primaryMapper.insertEquationData(primaryStateDot, yDot);
// Add contribution for secondary equations
for (final SecondaryComponent component : components) {
component.mapper.extractEquationData(y, component.state);
component.equation.computeDerivatives(t, primaryState, primaryStateDot,
component.state, component.stateDot);
component.mapper.insertEquationData(component.stateDot, yDot);
}
}
/** Add a set of secondary equations to be integrated along with the primary set.
* @param secondary secondary equations set
* @return index of the secondary equation in the expanded state
*/
public int addSecondaryEquations(final SecondaryEquations secondary) {
final int firstIndex;
if (components.isEmpty()) {
// lazy creation of the components list
components = new ArrayList<ExpandableStatefulODE.SecondaryComponent>();
firstIndex = primary.getDimension();
} else {
final SecondaryComponent last = components.get(components.size() - 1);
firstIndex = last.mapper.getFirstIndex() + last.mapper.getDimension();
}
components.add(new SecondaryComponent(secondary, firstIndex));
return components.size() - 1;
}
/** Get an equations mapper for the primary equations set.
* @return mapper for the primary set
* @see #getSecondaryMappers()
*/
public EquationsMapper getPrimaryMapper() {
return primaryMapper;
}
/** Get the equations mappers for the secondary equations sets.
* @return equations mappers for the secondary equations sets
* @see #getPrimaryMapper()
*/
public EquationsMapper[] getSecondaryMappers() {
final EquationsMapper[] mappers = new EquationsMapper[components.size()];
for (int i = 0; i < mappers.length; ++i) {
mappers[i] = components.get(i).mapper;
}
return mappers;
}
/** Set current time.
* @param time current time
*/
public void setTime(final double time) {
this.time = time;
}
/** Get current time.
* @return current time
*/
public double getTime() {
return time;
}
/** Set primary part of the current state.
* @param primaryState primary part of the current state
* @throws DimensionMismatchException if the dimension of the array does not
* match the primary set
*/
public void setPrimaryState(final double[] primaryState) throws DimensionMismatchException {
// safety checks
if (primaryState.length != this.primaryState.length) {
throw new DimensionMismatchException(primaryState.length, this.primaryState.length);
}
// set the data
System.arraycopy(primaryState, 0, this.primaryState, 0, primaryState.length);
}
/** Get primary part of the current state.
* @return primary part of the current state
*/
public double[] getPrimaryState() {
return primaryState.clone();
}
/** Get primary part of the current state derivative.
* @return primary part of the current state derivative
*/
public double[] getPrimaryStateDot() {
return primaryStateDot.clone();
}
/** Set secondary part of the current state.
* @param index index of the part to set as returned by {@link
* #addSecondaryEquations(SecondaryEquations)}
* @param secondaryState secondary part of the current state
* @throws DimensionMismatchException if the dimension of the partial state does not
* match the selected equations set dimension
*/
public void setSecondaryState(final int index, final double[] secondaryState)
throws DimensionMismatchException {
// get either the secondary state
double[] localArray = components.get(index).state;
// safety checks
if (secondaryState.length != localArray.length) {
throw new DimensionMismatchException(secondaryState.length, localArray.length);
}
// set the data
System.arraycopy(secondaryState, 0, localArray, 0, secondaryState.length);
}
/** Get secondary part of the current state.
* @param index index of the part to set as returned by {@link
* #addSecondaryEquations(SecondaryEquations)}
* @return secondary part of the current state
*/
public double[] getSecondaryState(final int index) {
return components.get(index).state.clone();
}
/** Get secondary part of the current state derivative.
* @param index index of the part to set as returned by {@link
* #addSecondaryEquations(SecondaryEquations)}
* @return secondary part of the current state derivative
*/
public double[] getSecondaryStateDot(final int index) {
return components.get(index).stateDot.clone();
}
/** Set the complete current state.
* @param completeState complete current state to copy data from
* @throws DimensionMismatchException if the dimension of the complete state does not
* match the complete equations sets dimension
*/
public void setCompleteState(final double[] completeState)
throws DimensionMismatchException {
// safety checks
if (completeState.length != getTotalDimension()) {
throw new DimensionMismatchException(completeState.length, getTotalDimension());
}
// set the data
primaryMapper.extractEquationData(completeState, primaryState);
for (final SecondaryComponent component : components) {
component.mapper.extractEquationData(completeState, component.state);
}
}
/** Get the complete current state.
* @return complete current state
* @throws DimensionMismatchException if the dimension of the complete state does not
* match the complete equations sets dimension
*/
public double[] getCompleteState() {
// allocate complete array
double[] completeState = new double[getTotalDimension()];
// set the data
primaryMapper.insertEquationData(primaryState, completeState);
for (final SecondaryComponent component : components) {
component.mapper.insertEquationData(component.state, completeState);
}
return completeState;
}
/** Components of the compound stateful ODE. */
private static class SecondaryComponent {
/** Secondary differential equation. */
private final SecondaryEquations equation;
/** Mapper between local and complete arrays. */
private final EquationsMapper mapper;
/** State. */
private final double[] state;
/** State derivative. */
private final double[] stateDot;
/** Simple constructor.
* @param equation secondary differential equation
* @param first index index to use for the first element in the complete arrays
*/
public SecondaryComponent(final SecondaryEquations equation, final int firstIndex) {
final int n = equation.getDimension();
this.equation = equation;
mapper = new EquationsMapper(firstIndex, n);
state = new double[n];
stateDot = new double[n];
}
}
}

66
src/main/java/org/apache/commons/math/ode/ExtendedFirstOrderDifferentialEquations.java
View File

@ -1,66 +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;
/** This interface represents a first order differential equations set
* with a main set of equations and an extension set.
*
* <p>
* This interface is a simple extension on the {@link
* FirstOrderDifferentialEquations} that allows to identify which part
* of a complete set of differential equations correspond to the main
* set and which part correspond to the extension set.
* </p>
* <p>
* One typical use case is the computation of Jacobians. The main
* set of equations correspond to the raw ode, and we add to this set
* another bunch of equations which represent the jacobians of the
* main set. In that case, we want the integrator to use <em>only</em>
* the main set to estimate the errors and hence the step sizes. It should
* <em>not</em> use the additional equations in this computation. If the
* complete ode implements this interface, the {@link FirstOrderIntegrator
* integrator} will be able to know where the main set ends and where the
* extended set begins.
* </p>
* <p>
* We consider that the main set always corresponds to the first equations
* and the extended set to the last equations.
* </p>
*
* @see FirstOrderDifferentialEquations
*
* @version $Id$
* @since 2.2
*/
public interface ExtendedFirstOrderDifferentialEquations extends FirstOrderDifferentialEquations {
/** Return the dimension of the main set of equations.
* <p>
* The main set of equations represent the first part of an ODE state.
* The error estimations and adaptive step size computation should be
* done on this first part only, not on the final part of the state
* which represent an extension set of equations which are considered
* secondary.
* </p>
* @return dimension of the main set of equations, must be lesser than or
* equal to the {@link #getDimension() total dimension}
*/
int getMainSetDimension();
}

601
src/main/java/org/apache/commons/math/ode/JacobianMatrices.java
View File

@ -25,25 +25,25 @@ import org.apache.commons.math.exception.MathIllegalArgumentException;
import org.apache.commons.math.exception.util.LocalizedFormats;
/**
* This class defines a set of {@link AdditionalEquations additional equations} to
* compute the jacobian matrices with respect to the initial state vector and, if
* any, to some parameters of the main ODE set.
* This class defines a set of {@link SecondaryEquations secondary equations} to
* compute the Jacobian matrices with respect to the initial state vector and, if
* any, to some parameters of the primary ODE set.
* <p>
* It is intended to be packed into an {@link ExpandableFirstOrderDifferentialEquations}
* in conjunction with a main set of ODE, which may be:
* It is intended to be packed into an {@link ExpandableStatefulODE}
* in conjunction with a primary set of ODE, which may be:
* <ul>
* <li>a {@link FirstOrderDifferentialEquations}</li>
* <li>a {@link MainStateJacobianProvider}</li>
* </ul>
* In order to compute jacobian matrices with respect to some parameters of the
* main ODE set, the following parameter jacobian providers may be set:
* In order to compute Jacobian matrices with respect to some parameters of the
* primary ODE set, the following parameter Jacobian providers may be set:
* <ul>
* <li>a {@link ParameterJacobianProvider}</li>
* <li>a {@link ParameterizedODE}</li>
* </ul>
* </p>
*
* @see ExpandableFirstOrderDifferentialEquations
* @see ExpandableStatefulODE
* @see FirstOrderDifferentialEquations
* @see MainStateJacobianProvider
* @see ParameterJacobianProvider
@ -52,288 +52,242 @@ import org.apache.commons.math.exception.util.LocalizedFormats;
* @version $Id$
* @since 3.0
*/
public class JacobianMatrices implements AdditionalEquations {
public class JacobianMatrices {
/** Expandable first order differential equation. */
private ExpandableFirstOrderDifferentialEquations efode;
private ExpandableStatefulODE efode;
/** FODE without exact main jacobian computation skill. */
private FirstOrderDifferentialEquations fode = null;
/** Index of the instance in the expandable set. */
private int index;
/** FODE with exact main jacobian computation skill. */
private MainStateJacobianProvider jode = null;
/** FODE with exact primary Jacobian computation skill. */
private MainStateJacobianProvider jode;
/** FODE without exact parameter jacobian computation skill. */
private ParameterizedODE pode = null;
/** FODE with exact parameter jacobian computation skill. */
private List<ParameterJacobianProvider> pjp = new ArrayList<ParameterJacobianProvider>();;
/** List of parameters selected for parameter jacobian computation. */
private List<ParameterConfiguration> selectedParameters = null;
/** FODE without exact parameter Jacobian computation skill. */
private ParameterizedODE pode;
/** Main state vector dimension. */
private int stateDim;
/** Selected parameters for parameter Jacobian computation. */
private ParameterConfiguration[] selectedParameters;
/** FODE with exact parameter Jacobian computation skill. */
private List<ParameterJacobianProvider> jacobianProviders;
/** Parameters dimension. */
private int paramDim = 0;
/** Current main state jacobian matrix in a row. */
private double[] mainJacobianInARow;
/** Current parameters jacobian matrices in a row. */
private double[] parameterJacobiansInARow = null;
/** Step used for finite difference computation of jacobian matrix
* w.r.t. main state vector. */
private double[] hY = null;
/** Boolean for fode consistency. */
private boolean dirtyMainState = false;
private int paramDim;
/** Boolean for selected parameters consistency. */
private boolean dirtyParameter = false;
private boolean dirtyParameter;
/** Simple constructor for an additional equations set computing jacobian matrices.
* <p>This additional equations set is added internally to the expandable
* first order differential equations set thanks to the
* {@link ExpandableFirstOrderDifferentialEquations#addAdditionalEquations(AdditionalEquations)}
* method.
* @param extended the expandable first order differential equations set
* @param jode the main first order differential equations set to extend
* @exception IllegalArgumentException if jode does not match the main set to be extended given by
* {@link ExpandableFirstOrderDifferentialEquations#getMainSet() extended.getMainSet()}
/** State and parameters Jacobian matrices in a row. */
private double[] matricesData;
/** Simple constructor for a secondary equations set computing Jacobian matrices.
* <p>
* Parameters must belong to the supported ones given by {@link
* Parameterizable#getParametersNames()}, so the primary set of differential
* equations must be {@link Parameterizable}.
* </p>
* <p>Note that each selection clears the previous selected parameters.</p>
*
* @param fode the primary first order differential equations set to extend
* @param hY step used for finite difference computation with respect to state vector
* @param parameters parameters to consider for Jacobian matrices processing
* (may be null if parameters Jacobians is not desired)
* @exception MathIllegalArgumentException if one parameter is not supported
* or there is a dimension mismatch with {@code hY}
*/
public JacobianMatrices(final ExpandableFirstOrderDifferentialEquations extended,
final MainStateJacobianProvider jode)
throws IllegalArgumentException {
public JacobianMatrices(final FirstOrderDifferentialEquations fode, final double[] hY,
final String... parameters)
throws MathIllegalArgumentException {
this(new MainStateJacobianWrapper(fode, hY), parameters);
}
checkCompatibility(extended, jode);
/** Simple constructor for a secondary equations set computing Jacobian matrices.
* <p>
* Parameters must belong to the supported ones given by {@link
* Parameterizable#getParametersNames()}, so the primary set of differential
* equations must be {@link Parameterizable}.
* </p>
* <p>Note that each selection clears the previous selected parameters.</p>
*
* @param jode the primary first order differential equations set to extend
* @param parameters parameters to consider for Jacobian matrices processing
* (may be null if parameters Jacobians is not desired)
* @exception MathIllegalArgumentException if one parameter is not supported
*/
public JacobianMatrices(final MainStateJacobianProvider jode,
final String... parameters)
throws MathIllegalArgumentException {
this.efode = null;
this.index = -1;
efode = extended;
stateDim = efode.getMainSetDimension();
mainJacobianInARow = new double[stateDim * stateDim];
this.jode = jode;
efode.addAdditionalEquations(this);
setInitialMainStateJacobian();
this.pode = null;
this.stateDim = jode.getDimension();
if (parameters == null) {
selectedParameters = null;
paramDim = 0;
} else {
this.selectedParameters = new ParameterConfiguration[parameters.length];
for (int i = 0; i < parameters.length; ++i) {
selectedParameters[i] = new ParameterConfiguration(parameters[i], Double.NaN);
}
paramDim = parameters.length;
}
this.dirtyParameter = false;
this.jacobianProviders = new ArrayList<ParameterJacobianProvider>();
// set the default initial state Jacobian to the identity
// and the default initial parameters Jacobian to the null matrix
matricesData = new double[(stateDim + paramDim) * stateDim];
for (int i = 0; i < stateDim; ++i) {
matricesData[i * (stateDim + 1)] = 1.0;
}
}
/** Simple constructor for an additional equations set computing jacobian matrices.
* <p>This additional equations set is added internally to the expandable
* first order differential equations set thanks to the
* {@link ExpandableFirstOrderDifferentialEquations#addAdditionalEquations(AdditionalEquations)}
* method.
* @param extended the expandable first order differential equations set
* @param fode the main first order differential equations set to extend
* @exception IllegalArgumentException if fode does not match the main set to be extended given by
* {@link ExpandableFirstOrderDifferentialEquations#getMainSet() extended.getMainSet()}
/** Register the variational equations for the Jacobians matrices to the expandable set.
* @exception MathIllegalArgumentException if the primary set of the expandable set does
* not match the one used to build the instance
* @see ExpandableStatefulODE#addSecondaryEquations(SecondaryEquations)
*/
public JacobianMatrices(final ExpandableFirstOrderDifferentialEquations extended,
final FirstOrderDifferentialEquations fode)
throws IllegalArgumentException {
public void registerVariationalEquations(final ExpandableStatefulODE expandable)
throws MathIllegalArgumentException {
checkCompatibility(extended, fode);
// safety checks
final FirstOrderDifferentialEquations ode = (jode instanceof MainStateJacobianWrapper) ?
((MainStateJacobianWrapper) jode).ode :
jode;
if (expandable.getPrimary() != ode) {
throw new MathIllegalArgumentException(LocalizedFormats.UNMATCHED_ODE_IN_EXPANDED_SET);
}
efode = expandable;
index = efode.addSecondaryEquations(new JacobiansSecondaryEquations());
efode.setSecondaryState(index, matricesData);
efode = extended;
stateDim = efode.getMainSetDimension();
mainJacobianInARow = new double[stateDim * stateDim];
this.fode = fode;
dirtyMainState = true;
efode.addAdditionalEquations(this);
setInitialMainStateJacobian();
}
/** Add a parameter jacobian provider.
* @param pjp the parameter jacobian provider to compute exactly the parameter jacobian matrix
/** Add a parameter Jacobian provider.
* @param provider the parameter Jacobian provider to compute exactly the parameter Jacobian matrix
*/
public void setParameterJacobianProvider(final ParameterJacobianProvider pjp) {
this.pjp.add(pjp);
public void addParameterJacobianProvider(final ParameterJacobianProvider provider) {
jacobianProviders.add(provider);
}
/** Add a parameter jacobian provider.
* @param pjp the parameterized ODE to compute by finite difference the parameter jacobian matrix
/** Add a parameter Jacobian provider.
* @param pode the parameterized ODE to compute the parameter Jacobian matrix using finite differences
*/
public void setParameterizedODE(final ParameterizedODE pode) {
this.pode = pode;
dirtyParameter = true;
}
/** Select the parameters to consider for jacobian matrices processing.
* <p>
* Parameters must belong to the supported ones given by {@link
* Parameterizable#getParametersNames()}, so the main set of differential
* equations must be {@link Parameterizable}.
* </p>
* <p>Note that each selection clears the previous selected parameters.</p>
*
* @param parameters parameters to consider for jacobian matrices processing
* @exception IllegalArgumentException if one parameter is not supported
*/
public void selectParameters(final String... parameters) throws IllegalArgumentException {
selectedParameters = new ArrayList<ParameterConfiguration>();
for (String param : parameters) {
selectedParameters.add(new ParameterConfiguration(param, Double.NaN));
}
paramDim = parameters.length;
parameterJacobiansInARow = new double[paramDim * stateDim];
setInitialParameterJacobians();
}
/** Set the step associated to a parameter in order to compute by finite
* difference the jacobian matrix.
* difference the Jacobian matrix.
* <p>
* Needed if and only if the main ODE set is a {@link ParameterizedODE}
* and the parameter has been {@link #selectParameters(String ...) selected}
* Needed if and only if the primary ODE set is a {@link ParameterizedODE}.
* </p>
* <p>
* For pval, a non zero value of the parameter, pval * Math.sqrt(MathUtils.EPSILON)
* is a reasonable value for such a step.
* Given a non zero parameter value pval for the parameter, a reasonable value
* for such a step is {@code pval * FastMath.sqrt(MathUtils.EPSILON)}.
* </p>
* <p>
* A zero value for such a step doesn't enable to compute the parameter jacobian matrix.
* A zero value for such a step doesn't enable to compute the parameter Jacobian matrix.
* </p>
* @param parameter parameter to consider for jacobian processing
* @param hP step for jacobian finite difference computation w.r.t. the specified parameter
* @param parameter parameter to consider for Jacobian processing
* @param hP step for Jacobian finite difference computation w.r.t. the specified parameter
* @see ParameterizedODE
* @exception IllegalArgumentException if the parameter is not supported
*/
public void setParameterStep(final String parameter, final double hP) {
boolean found = false;
for (ParameterConfiguration param: selectedParameters) {
if (parameter.equals(param.getParameterName())) {
param.setHP(hP);
found = true;
dirtyParameter = true;
break;
return;
}
}
if (!found) {
throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER,
parameter);
}
throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER, parameter);
}
/** Set the steps in order to compute by finite difference the jacobian
* matrix with respect to main state.
/** Set the initial value of the Jacobian matrix with respect to state.
* <p>
* Needed if and only if the main set is a {@link FirstOrderDifferentialEquations}.
* If this method is not called, the initial value of the Jacobian
* matrix with respect to state is set to identity.
* </p>
* <p>
* Zero values for those steps don't enable to compute the main state jacobian matrix.
* </p>
* @param hY step used for finite difference computation with respect to state vector
* @exception IllegalArgumentException if the hY has not the dimension of the main state
* given by {@link ExpandableFirstOrderDifferentialEquations#getMainSetDimension()}
*/
public void setMainStateSteps(final double[] hY) {
if (fode != null) {
// Check dimension
checkDimension(stateDim, hY);
this.hY = hY.clone();
dirtyMainState = true;
}
}
/** Set the initial value of the jacobian matrix with respect to state.
* @param dYdY0 initial jacobian matrix w.r.t. state
* @param dYdY0 initial Jacobian matrix w.r.t. state
* @exception IllegalArgumentException if matrix dimensions are incorrect
*/
public void setInitialMainStateJacobian(final double[][] dYdY0) {
public void setInitialMainStateJacobian(final double[][] dYdY0)
throws MathIllegalArgumentException {
// Check dimensions
checkDimension(stateDim, dYdY0);
checkDimension(stateDim, dYdY0[0]);
// store the matrix in row major order as a single dimension array
int index = 0;
int i = 0;
for (final double[] row : dYdY0) {
System.arraycopy(row, 0, mainJacobianInARow, index, stateDim);
index += stateDim;
System.arraycopy(row, 0, matricesData, i, stateDim);
i += stateDim;
}
// set initial additional state value in expandable fode
efode.setInitialAdditionalState(mainJacobianInARow, this);
if (efode != null) {
efode.setSecondaryState(index, matricesData);
}
}
/** Set the initial value of the jacobian matrix with respect to one parameter.
* <p>The parameter must be {@link #selectParameters(String...) selected}.</p>
/** Set the initial value of a column of the Jacobian matrix with respect to one parameter.
* <p>
* If this method is not called for some parameter, the initial value of
* the column of the Jacobian matrix with respect to this parameter is set to zero.
* </p>
* @param pName parameter name
* @param dYdP initial jacobian matrix w.r.t. the parameter
* @exception IllegalArgumentException if matrix dimensions are incorrect
* @param dYdP initial Jacobian column vector with respect to the parameter
* @exception MathIllegalArgumentException if a parameter is not supported
*/
public void setInitialParameterJacobian(final String pName, final double[] dYdP) {
public void setInitialParameterJacobian(final String pName, final double[] dYdP)
throws MathIllegalArgumentException {
// Check dimensions
checkDimension(stateDim, dYdP);
// store the matrix in a global single dimension array
boolean found = false;
int index = 0;
// store the column in a global single dimension array
int i = stateDim * stateDim;
for (ParameterConfiguration param: selectedParameters) {
if (pName.equals(param.getParameterName())) {
System.arraycopy(dYdP, 0, parameterJacobiansInARow, index, stateDim);
double[] p = new double[this.getDimension()];
index = stateDim * stateDim;
System.arraycopy(mainJacobianInARow, 0, p, 0, index);
System.arraycopy(parameterJacobiansInARow, 0, p, index, stateDim * paramDim);
// set initial additional state value in expandable fode
efode.setInitialAdditionalState(p, this);
found = true;
break;
System.arraycopy(dYdP, 0, matricesData, i, stateDim);
if (efode != null) {
efode.setSecondaryState(index, matricesData);
}
return;
}
index += stateDim;
}
if (! found) {
throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER,
pName);
i += stateDim;
}
throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER, pName);
}
/** Set the default initial value of the jacobian matrix with respect to state.
* <p>dYdY0 is set to the identity matrix.</p>
*/
public void setInitialMainStateJacobian() {
final double[][] dYdY0 = new double[stateDim][stateDim];
for (int i = 0; i < stateDim; ++i) {
dYdY0[i][i] = 1.0;
}
setInitialMainStateJacobian(dYdY0);
}
/** Set the default initial value of the jacobian matrix with respect to one parameter.
* <p>The parameter must be {@link #selectParameters(String...) selected}.</p>
* <p>dYdP is set to the null matrix.</p>
* @param pName parameter name
*/
public void setInitialParameterJacobian(final String pName) {
setInitialParameterJacobian(pName, new double[stateDim]);
}
/** Set the default initial values of jacobian matrices with respect to all parameters.
*/
public void setInitialParameterJacobians() {
for (ParameterConfiguration param: selectedParameters) {
setInitialParameterJacobian(param.getParameterName());
}
}
/** Set default initial values for jacobian matrices.
* <p>dYdY0 is set to the identity matrix and all dYdP are set to zero.</p>
*/
public void setInitialJacobians() {
setInitialMainStateJacobian();
setInitialParameterJacobians();
}
/** Get the current value of the jacobian matrix with respect to state.
* @param dYdY0 current jacobian matrix with respect to state.
/** Get the current value of the Jacobian matrix with respect to state.
* @param dYdY0 current Jacobian matrix with respect to state.
*/
public void getCurrentMainSetJacobian(final double[][] dYdY0) {
// get current state for this set of equations from the expandable fode
double[] p = efode.getCurrentAdditionalState(this);
double[] p = efode.getSecondaryState(index);
int index = 0;
for (int i = 0; i < stateDim; i++) {
@ -343,14 +297,14 @@ public class JacobianMatrices implements AdditionalEquations {
}
/** Get the current value of the jacobian matrix with respect to one parameter.
* @param pName name of the parameter for the computed jacobian matrix
* @param dYdP current jacobian matrix with respect to the named parameter
/** Get the current value of the Jacobian matrix with respect to one parameter.
* @param pName name of the parameter for the computed Jacobian matrix
* @param dYdP current Jacobian matrix with respect to the named parameter
*/
public void getCurrentParameterJacobian(String pName, final double[] dYdP) {
// get current state for this set of equations from the expandable fode
double[] p = efode.getCurrentAdditionalState(this);
double[] p = efode.getSecondaryState(index);
int index = stateDim * stateDim;
for (ParameterConfiguration param: selectedParameters) {
@ -363,97 +317,6 @@ public class JacobianMatrices implements AdditionalEquations {
}
/** {@inheritDoc} */
public int getDimension() {
return stateDim * (stateDim + paramDim);
}
/** {@inheritDoc} */
public void computeDerivatives(final double t, final double[] y, final double[] yDot,
final double[] z, final double[] zDot) {
// Lazy initialization
if (dirtyMainState) {
jode = new MainStateJacobianWrapper(fode, hY);
dirtyMainState = false;
}
if (dirtyParameter && (paramDim != 0)) {
pjp.add(new ParameterJacobianWrapper(jode, pode, selectedParameters));
dirtyParameter = false;
}
// variational equations:
// from d[dy/dt]/dy0 and d[dy/dt]/dp to d[dy/dy0]/dt and d[dy/dp]/dt
// compute jacobian matrix with respect to main state
double[][] dFdY = new double[stateDim][stateDim];
jode.computeMainStateJacobian(t, y, yDot, dFdY);
// Dispatch jacobian matrix in the compound additional state vector
for (int i = 0; i < stateDim; ++i) {
final double[] dFdYi = dFdY[i];
for (int j = 0; j < stateDim; ++j) {
double s = 0;
final int startIndex = j;
int zIndex = startIndex;
for (int l = 0; l < stateDim; ++l) {
s += dFdYi[l] * z[zIndex];
zIndex += stateDim;
}
zDot[startIndex + i * stateDim] = s;
}
}
if (paramDim != 0) {
// compute jacobian matrices with respect to parameters
double[] dFdP = new double[stateDim];
int startIndex = stateDim * stateDim;
for (ParameterConfiguration param: selectedParameters) {
boolean found = false;
for (ParameterJacobianProvider provider: pjp) {
if (provider.isSupported(param.getParameterName())) {
try {
provider.computeParameterJacobian(t, y, yDot, param.getParameterName(), dFdP);
for (int i = 0; i < stateDim; ++i) {
final double[] dFdYi = dFdY[i];
int zIndex = startIndex;
double s = dFdP[i];
for (int l = 0; l < stateDim; ++l) {
s += dFdYi[l] * z[zIndex];
zIndex++;
}
zDot[startIndex + i] = s;
}
startIndex += stateDim;
found = true;
break;
} catch (IllegalArgumentException iae) {
}
}
}
if (! found) {
throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER,
param);
}
}
}
}
/** Check compatibility between the main set in the expandable ode and an ordinary ode.
* @param expended expandable ode containing a main set
* @param ode single ode to check
* @throws MathIllegalArgumentException if single ode doesn't match the main ode set in the extended ode
*/
private void checkCompatibility(final ExpandableFirstOrderDifferentialEquations extended,
final FirstOrderDifferentialEquations ode)
throws MathIllegalArgumentException {
if (!(ode == extended.getMainSet())) {
throw new MathIllegalArgumentException(LocalizedFormats.UNMATCHED_ODE_IN_EXTENDED_SET);
}
}
/** Check array dimensions.
* @param expected expected dimension
* @param array (may be null if expected is 0)
@ -467,5 +330,139 @@ public class JacobianMatrices implements AdditionalEquations {
}
}
/** Local implementation of secondary equations.
* <p>
* This class is an inner class to ensure proper scheduling of calls
* by forcing the use of {@link JacobianMatrices#registerVariationalEquations(ExpandableStatefulODE)}.
* </p>
*/
private class JacobiansSecondaryEquations implements SecondaryEquations {
/** {@inheritDoc} */
public int getDimension() {
return stateDim * (stateDim + paramDim);
}
/** {@inheritDoc} */
public void computeDerivatives(final double t, final double[] y, final double[] yDot,
final double[] z, final double[] zDot) {
// Lazy initialization
if (dirtyParameter && (paramDim != 0)) {
jacobianProviders.add(new ParameterJacobianWrapper(jode, pode, selectedParameters));
dirtyParameter = false;
}
// variational equations:
// from d[dy/dt]/dy0 and d[dy/dt]/dp to d[dy/dy0]/dt and d[dy/dp]/dt
// compute Jacobian matrix with respect to primary state
double[][] dFdY = new double[stateDim][stateDim];
jode.computeMainStateJacobian(t, y, yDot, dFdY);
// Dispatch Jacobian matrix in the compound secondary state vector
for (int i = 0; i < stateDim; ++i) {
final double[] dFdYi = dFdY[i];
for (int j = 0; j < stateDim; ++j) {
double s = 0;
final int startIndex = j;
int zIndex = startIndex;
for (int l = 0; l < stateDim; ++l) {
s += dFdYi[l] * z[zIndex];
zIndex += stateDim;
}
zDot[startIndex + i * stateDim] = s;
}
}
if (paramDim != 0) {
// compute Jacobian matrices with respect to parameters
double[] dFdP = new double[stateDim];
int startIndex = stateDim * stateDim;
for (ParameterConfiguration param: selectedParameters) {
boolean found = false;
for (ParameterJacobianProvider provider: jacobianProviders) {
if (provider.isSupported(param.getParameterName())) {
try {
provider.computeParameterJacobian(t, y, yDot, param.getParameterName(), dFdP);
for (int i = 0; i < stateDim; ++i) {
final double[] dFdYi = dFdY[i];
int zIndex = startIndex;
double s = dFdP[i];
for (int l = 0; l < stateDim; ++l) {
s += dFdYi[l] * z[zIndex];
zIndex++;
}
zDot[startIndex + i] = s;
}
startIndex += stateDim;
found = true;
break;
} catch (IllegalArgumentException iae) {
}
}
}
if (! found) {
throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER,
param);
}
}
}
}
}
/** Wrapper class to compute jacobian matrices by finite differences for ODE
* which do not compute them by themselves.
*/
private static class MainStateJacobianWrapper implements MainStateJacobianProvider {
/** Raw ODE without jacobians computation skill to be wrapped into a MainStateJacobianProvider. */
private final FirstOrderDifferentialEquations ode;
/** Steps for finite difference computation of the jacobian df/dy w.r.t. state. */
private final double[] hY;
/** Wrap a {@link FirstOrderDifferentialEquations} into a {@link MainStateJacobianProvider}.
* @param ode original ODE problem, without jacobians computation skill
* @param hY step sizes to compute the jacobian df/dy
* @see JacobianMatrices#setMainStateSteps(double[])
*/
public MainStateJacobianWrapper(final FirstOrderDifferentialEquations ode,
final double[] hY) {
this.ode = ode;
this.hY = hY.clone();
}
/** {@inheritDoc} */
public int getDimension() {
return ode.getDimension();
}
/** {@inheritDoc} */
public void computeDerivatives(double t, double[] y, double[] yDot) {
ode.computeDerivatives(t, y, yDot);
}
/** {@inheritDoc} */
public void computeMainStateJacobian(double t, double[] y, double[] yDot,
double[][] dFdY) {
final int n = ode.getDimension();
final double[] tmpDot = new double[n];
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;
}
}
}
}

72
src/main/java/org/apache/commons/math/ode/MainStateJacobianWrapper.java
View File

@ -1,72 +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;
/** Wrapper class to compute jacobian matrices by finite differences for ODE
* which do not compute them by themselves.
*
* @version $Id$
* @since 3.0
*/
class MainStateJacobianWrapper implements MainStateJacobianProvider {
/** Raw ODE without jacobians computation skill to be wrapped into a MainStateJacobianProvider. */
private final FirstOrderDifferentialEquations ode;
/** Steps for finite difference computation of the jacobian df/dy w.r.t. state. */
private final double[] hY;
/** Wrap a {@link FirstOrderDifferentialEquations} into a {@link MainStateJacobianProvider}.
* @param ode original ODE problem, without jacobians computation skill
* @param hY step sizes to compute the jacobian df/dy
* @see JacobianMatrices#setMainStateSteps(double[])
*/
public MainStateJacobianWrapper(final FirstOrderDifferentialEquations ode,
final double[] hY) {
this.ode = ode;
this.hY = hY.clone();
}
/** {@inheritDoc} */
public int getDimension() {
return ode.getDimension();
}
/** {@inheritDoc} */
public void computeDerivatives(double t, double[] y, double[] yDot) {
ode.computeDerivatives(t, y, yDot);
}
/** {@inheritDoc} */
public void computeMainStateJacobian(double t, double[] y, double[] yDot,
double[][] dFdY) {
final int n = ode.getDimension();
final double[] tmpDot = new double[n];
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;
}
}
}

7
src/main/java/org/apache/commons/math/ode/ParameterJacobianProvider.java
View File

@ -18,7 +18,7 @@ package org.apache.commons.math.ode;
import org.apache.commons.math.exception.MathIllegalArgumentException;
/** Interface to compute exactly jacobian matrix for some parameter
/** Interface to compute exactly Jacobian matrix for some parameter
* when computing {@link JacobianMatrices partial derivatives equations}.
*
* @version $Id$
@ -26,18 +26,19 @@ import org.apache.commons.math.exception.MathIllegalArgumentException;
*/
public interface ParameterJacobianProvider extends Parameterizable {
/** Compute the jacobian matrix of ODE with respect to one parameter.
/** Compute the Jacobian matrix of ODE with respect to one parameter.
* <p>The parameter must be one given by {@link #getParametersNames()}.</p>
* @param t current value of the independent <I>time</I> variable
* @param y array containing the current value of the main state vector
* @param yDot array containing the current value of the time derivative
* of the main state vector
* @param paramName name of the parameter to consider
* @param dFdP placeholder array where to put the jacobian matrix of the
* @param dFdP placeholder array where to put the Jacobian matrix of the
* ODE with respect to the parameter
* @throws MathIllegalArgumentException if the parameter is not supported
*/
void computeParameterJacobian(double t, double[] y, double[] yDot,
String paramName, double[] dFdP)
throws MathIllegalArgumentException;
}

12
src/main/java/org/apache/commons/math/ode/ParameterJacobianWrapper.java
View File

@ -20,7 +20,7 @@ import java.util.Collection;
import java.util.HashMap;
import java.util.Map;
/** Wrapper class to compute jacobian matrices by finite differences for ODE
/** Wrapper class to compute Jacobian matrices by finite differences for ODE
* which do not compute them by themselves.
*
* @version $Id$
@ -31,21 +31,21 @@ class ParameterJacobianWrapper implements ParameterJacobianProvider {
/** Main ODE set. */
private final FirstOrderDifferentialEquations fode;
/** Raw ODE without jacobian computation skill to be wrapped into a ParameterJacobianProvider. */
/** Raw ODE without Jacobian computation skill to be wrapped into a ParameterJacobianProvider. */
private final ParameterizedODE pode;
/** Steps for finite difference computation of the jacobian df/dp w.r.t. parameters. */
/** Steps for finite difference computation of the Jacobian df/dp w.r.t. parameters. */
private final Map<String, Double> hParam;
/** Wrap a {@link ParameterizedODE} into a {@link ParameterJacobianProvider}.
* @param fode main first order differential equations set
* @param pode additional problem, without parametre jacobian computation skill
* @param paramsAndSteps parameters and steps to compute the jacobians df/dp
* @param pode secondary problem, without parameter Jacobian computation skill
* @param paramsAndSteps parameters and steps to compute the Jacobians df/dp
* @see JacobianMatrices#setParameterStep(String, double)
*/
public ParameterJacobianWrapper(final FirstOrderDifferentialEquations fode,
final ParameterizedODE pode,
final Collection<ParameterConfiguration> paramsAndSteps) {
final ParameterConfiguration[] paramsAndSteps) {
this.fode = fode;
this.pode = pode;
this.hParam = new HashMap<String, Double>();

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

@ -16,7 +16,7 @@
*/
package org.apache.commons.math.ode;
/** Interface to compute by finite difference jacobian matrix for some parameter
/** Interface to compute by finite difference Jacobian matrix for some parameter
* when computing {@link JacobianMatrices partial derivatives equations}.
*
* @version $Id$

29
src/main/java/org/apache/commons/math/ode/AdditionalEquations.java → src/main/java/org/apache/commons/math/ode/SecondaryEquations.java
View File

@ -18,38 +18,39 @@
package org.apache.commons.math.ode;
/**
* This interface allows users to add their own differential equations to a main
* This interface allows users to add secondary differential equations to a primary
* set of differential equations.
* <p>
* In some cases users may need to integrate some problem-specific equations along
* with a main set of differential equations. One example is optimal control where
* with a primary set of differential equations. One example is optimal control where
* adjoined parameters linked to the minimized hamiltonian must be integrated.
* </p>
* <p>
* This interface allows users to add such equations to a main set of {@link
* This interface allows users to add such equations to a primary set of {@link
* FirstOrderDifferentialEquations first order differential equations}
* thanks to the {@link
* ExpandableFirstOrderDifferentialEquations#addAdditionalEquations(AdditionalEquations)}
* ExpandableStatefulODE#addSecondaryEquations(SecondaryEquations)}
* method.
* </p>
* @see ExpandableFirstOrderDifferentialEquations
* @see ExpandableStatefulODE
* @version $Id$
* @since 3.0
*/
public interface AdditionalEquations {
public interface SecondaryEquations {
/** Get the dimension of the additional state parameters.
* @return dimension of the additional state parameters
/** Get the dimension of the secondary state parameters.
* @return dimension of the secondary state parameters
*/
int getDimension();
/** Compute the derivatives related to the additional state parameters.
/** Compute the derivatives related to the secondary state parameters.
* @param t current value of the independent <I>time</I> variable
* @param y array containing the current value of the main state vector
* @param yDot array containing the derivative of the main state vector
* @param z array containing the current value of the additional state vector
* @param zDot placeholder array where to put the derivative of the additional state vector
* @param primary array containing the current value of the primary state vector
* @param primaryDot array containing the derivative of the primary state vector
* @param secondary array containing the current value of the secondary state vector
* @param secondaryDot placeholder array where to put the derivative of the secondary state vector
*/
void computeDerivatives(double t, double[] y, double[] yDot, double[] z, double[] zDot);
void computeDerivatives(double t, double[] primary, double[] primaryDot,
double[] secondary, double[] secondaryDot);
}

36
src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
View File

@ -20,7 +20,7 @@ package org.apache.commons.math.ode.nonstiff;
import org.apache.commons.math.exception.MathIllegalArgumentException;
import org.apache.commons.math.exception.MathIllegalStateException;
import org.apache.commons.math.linear.Array2DRowRealMatrix;
import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations;
import org.apache.commons.math.ode.ExpandableStatefulODE;
import org.apache.commons.math.ode.sampling.NordsieckStepInterpolator;
import org.apache.commons.math.ode.sampling.StepHandler;
import org.apache.commons.math.util.FastMath;
@ -189,31 +189,23 @@ public class AdamsBashforthIntegrator extends AdamsIntegrator {
/** {@inheritDoc} */
@Override
public double integrate(final ExpandableFirstOrderDifferentialEquations equations,
final double t0, final double[] z0,
final double t, final double[] z)
public void integrate(final ExpandableStatefulODE equations, final double t)
throws MathIllegalStateException, MathIllegalArgumentException {
sanityChecks(equations, t0, z0, t, z);
sanityChecks(equations, t);
setEquations(equations);
resetEvaluations();
final boolean forward = t > t0;
final boolean forward = t > equations.getTime();
// initialize working arrays
final int totalDim = equations.getDimension();
final int mainDim = equations.getMainSetDimension();
final double[] y0 = new double[totalDim];
final double[] y = new double[totalDim];
System.arraycopy(z0, 0, y0, 0, mainDim);
System.arraycopy(equations.getCurrentAdditionalStates(), 0, y0, mainDim, totalDim - mainDim);
if (y != y0) {
System.arraycopy(y0, 0, y, 0, totalDim);
}
final double[] yDot = new double[totalDim];
final double[] y0 = equations.getCompleteState();
final double[] y = y0.clone();
final double[] yDot = new double[y.length];
// set up an interpolator sharing the integrator arrays
final NordsieckStepInterpolator interpolator = new NordsieckStepInterpolator();
interpolator.reinitialize(y, forward);
interpolator.reinitialize(y, forward,
equations.getPrimaryMapper(), equations.getSecondaryMappers());
// set up integration control objects
for (StepHandler handler : stepHandlers) {
@ -222,7 +214,7 @@ public class AdamsBashforthIntegrator extends AdamsIntegrator {
setStateInitialized(false);
// compute the initial Nordsieck vector using the configured starter integrator
start(t0, y, t);
start(equations.getTime(), y, t);
interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
interpolator.storeTime(stepStart);
final int lastRow = nordsieck.getRowDimension() - 1;
@ -317,13 +309,11 @@ public class AdamsBashforthIntegrator extends AdamsIntegrator {
} while (!isLastStep);
// dispatch result between main and additional states
System.arraycopy(y, 0, z, 0, z.length);
equations.setCurrentAdditionalState(y);
// dispatch results
equations.setTime(stepStart);
equations.setCompleteState(y);
final double stopTime = stepStart;
resetInternalState();
return stopTime;
}

6
src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java
View File

@ -20,7 +20,7 @@ package org.apache.commons.math.ode.nonstiff;
import org.apache.commons.math.exception.MathIllegalArgumentException;
import org.apache.commons.math.exception.MathIllegalStateException;
import org.apache.commons.math.linear.Array2DRowRealMatrix;
import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations;
import org.apache.commons.math.ode.ExpandableStatefulODE;
import org.apache.commons.math.ode.MultistepIntegrator;
@ -86,9 +86,7 @@ public abstract class AdamsIntegrator extends MultistepIntegrator {
/** {@inheritDoc} */
@Override
public abstract double integrate(final ExpandableFirstOrderDifferentialEquations equations,
final double t0, final double[] y0,
final double t, final double[] y)
public abstract void integrate(final ExpandableStatefulODE equations, final double t)
throws MathIllegalStateException, MathIllegalArgumentException;
/** {@inheritDoc} */

45
src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
View File

@ -23,7 +23,7 @@ import org.apache.commons.math.exception.MathIllegalArgumentException;
import org.apache.commons.math.exception.MathIllegalStateException;
import org.apache.commons.math.linear.Array2DRowRealMatrix;
import org.apache.commons.math.linear.RealMatrixPreservingVisitor;
import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations;
import org.apache.commons.math.ode.ExpandableStatefulODE;
import org.apache.commons.math.ode.sampling.NordsieckStepInterpolator;
import org.apache.commons.math.ode.sampling.StepHandler;
import org.apache.commons.math.util.FastMath;
@ -206,34 +206,26 @@ public class AdamsMoultonIntegrator extends AdamsIntegrator {
/** {@inheritDoc} */
@Override
public double integrate(final ExpandableFirstOrderDifferentialEquations equations,
final double t0, final double[] z0,
final double t, final double[] z)
public void integrate(final ExpandableStatefulODE equations,final double t)
throws MathIllegalStateException, MathIllegalArgumentException {
sanityChecks(equations, t0, z0, t, z);
sanityChecks(equations, t);
setEquations(equations);
resetEvaluations();
final boolean forward = t > t0;
final boolean forward = t > equations.getTime();
// initialize working arrays
final int totalDim = equations.getDimension();
final int mainDim = equations.getMainSetDimension();
final double[] y0 = new double[totalDim];
final double[] y = new double[totalDim];
System.arraycopy(z0, 0, y0, 0, mainDim);
System.arraycopy(equations.getCurrentAdditionalStates(), 0, y0, mainDim, totalDim - mainDim);
if (y != y0) {
System.arraycopy(y0, 0, y, 0, totalDim);
}
final double[] yDot = new double[totalDim];
final double[] yTmp = new double[totalDim];
final double[] predictedScaled = new double[totalDim];
final double[] y0 = equations.getCompleteState();
final double[] y = y0.clone();
final double[] yDot = new double[y.length];
final double[] yTmp = new double[y.length];
final double[] predictedScaled = new double[y.length];
Array2DRowRealMatrix nordsieckTmp = null;
// set up two interpolators sharing the integrator arrays
final NordsieckStepInterpolator interpolator = new NordsieckStepInterpolator();
interpolator.reinitialize(y, forward);
interpolator.reinitialize(y, forward,
equations.getPrimaryMapper(), equations.getSecondaryMappers());
// set up integration control objects
for (StepHandler handler : stepHandlers) {
@ -242,7 +234,7 @@ public class AdamsMoultonIntegrator extends AdamsIntegrator {
setStateInitialized(false);
// compute the initial Nordsieck vector using the configured starter integrator
start(t0, y, t);
start(equations.getTime(), y, t);
interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
interpolator.storeTime(stepStart);
@ -295,7 +287,7 @@ public class AdamsMoultonIntegrator extends AdamsIntegrator {
updateHighOrderDerivativesPhase2(predictedScaled, correctedScaled, nordsieckTmp);
// discrete events handling
System.arraycopy(yTmp, 0, y, 0, totalDim);
System.arraycopy(yTmp, 0, y, 0, y.length);
interpolator.reinitialize(stepEnd, stepSize, correctedScaled, nordsieckTmp);
interpolator.storeTime(stepStart);
interpolator.shift();
@ -335,14 +327,11 @@ public class AdamsMoultonIntegrator extends AdamsIntegrator {
} while (!isLastStep);
// dispatch result between main and additional states
System.arraycopy(y, 0, z, 0, z.length);
equations.setCurrentAdditionalState(y);
// dispatch results
equations.setTime(stepStart);
equations.setCompleteState(y);
final double stopTime = stepStart;
stepStart = Double.NaN;
stepSize = Double.NaN;
return stopTime;
resetInternalState();
}

28
src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java
View File

@ -23,7 +23,7 @@ import org.apache.commons.math.exception.MathIllegalStateException;
import org.apache.commons.math.exception.NumberIsTooSmallException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.ode.AbstractIntegrator;
import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations;
import org.apache.commons.math.ode.ExpandableStatefulODE;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.util.FastMath;
@ -43,9 +43,9 @@ import org.apache.commons.math.util.FastMath;
* relTol which will be used for all components.
* </p>
* <p>
* If the Ordinary Differential Equations is an {@link ExpandableFirstOrderDifferentialEquations
* If the Ordinary Differential Equations is an {@link ExpandableStatefulODE
* extended ODE} rather than a {@link FirstOrderDifferentialEquations basic ODE}, then
* <em>only</em> the {@link ExpandableFirstOrderDifferentialEquations#getMainSet() main part}
* <em>only</em> the {@link ExpandableStatefulODE#getMainSet() main part}
* of the state vector is used for stepsize control, not the complete state vector.
* </p>
*
@ -213,24 +213,14 @@ public abstract class AdaptiveStepsizeIntegrator
}
}
/** Perform some sanity checks on the integration parameters.
* @param equations differential equations set
* @param t0 start time
* @param y0 state vector at t0
* @param t target time for the integration
* @param y placeholder where to put the state vector
* @exception DimensionMismatchException if some inconsistency is detected
* @exception NumberIsTooSmallException if integration span is too small
*/
/** {@inheritDoc} */
@Override
protected void sanityChecks(final ExpandableFirstOrderDifferentialEquations equations,
final double t0, final double[] y0,
final double t, final double[] y)
protected void sanityChecks(final ExpandableStatefulODE equations, final double t)
throws DimensionMismatchException, NumberIsTooSmallException {
super.sanityChecks(equations, t0, y0, t, y);
super.sanityChecks(equations, t);
mainSetDimension = equations.getMainSetDimension();
mainSetDimension = equations.getPrimaryMapper().getDimension();
if ((vecAbsoluteTolerance != null) && (vecAbsoluteTolerance.length != mainSetDimension)) {
throw new DimensionMismatchException(mainSetDimension, vecAbsoluteTolerance.length);
@ -349,9 +339,7 @@ public abstract class AdaptiveStepsizeIntegrator
}
/** {@inheritDoc} */
public abstract double integrate (ExpandableFirstOrderDifferentialEquations equations,
double t0, double[] y0,
double t, double[] y)
public abstract void integrate (ExpandableStatefulODE equations, double t)
throws MathIllegalStateException, MathIllegalArgumentException;
/** {@inheritDoc} */

7
src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java
View File

@ -18,6 +18,7 @@
package org.apache.commons.math.ode.nonstiff;
import org.apache.commons.math.ode.AbstractIntegrator;
import org.apache.commons.math.ode.EquationsMapper;
import org.apache.commons.math.ode.sampling.StepInterpolator;
/**
@ -145,8 +146,10 @@ class DormandPrince54StepInterpolator
/** {@inheritDoc} */
@Override
public void reinitialize(final AbstractIntegrator integrator,
final double[] y, final double[][] yDotK, final boolean forward) {
super.reinitialize(integrator, y, yDotK, forward);
final double[] y, final double[][] yDotK, final boolean forward,
final EquationsMapper primaryMapper,
final EquationsMapper[] secondaryMappers) {
super.reinitialize(integrator, y, yDotK, forward, primaryMapper, secondaryMappers);
v1 = null;
v2 = null;
v3 = null;

9
src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java
View File

@ -22,6 +22,7 @@ import java.io.ObjectInput;
import java.io.ObjectOutput;
import org.apache.commons.math.ode.AbstractIntegrator;
import org.apache.commons.math.ode.EquationsMapper;
import org.apache.commons.math.ode.sampling.StepInterpolator;
/**
@ -280,9 +281,11 @@ class DormandPrince853StepInterpolator
/** {@inheritDoc} */
@Override
public void reinitialize(final AbstractIntegrator integrator,
final double[] y, final double[][] yDotK, final boolean forward) {
final double[] y, final double[][] yDotK, final boolean forward,
final EquationsMapper primaryMapper,
final EquationsMapper[] secondaryMappers) {
super.reinitialize(integrator, y, yDotK, forward);
super.reinitialize(integrator, y, yDotK, forward, primaryMapper, secondaryMappers);
final int dimension = currentState.length;
@ -454,7 +457,7 @@ class DormandPrince853StepInterpolator
/** {@inheritDoc} */
@Override
public void readExternal(final ObjectInput in)
throws IOException {
throws IOException, ClassNotFoundException {
// read the local attributes
yDotKLast = new double[3][];

42
src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
View File

@ -19,7 +19,7 @@ package org.apache.commons.math.ode.nonstiff;
import org.apache.commons.math.exception.MathIllegalArgumentException;
import org.apache.commons.math.exception.MathIllegalStateException;
import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations;
import org.apache.commons.math.ode.ExpandableStatefulODE;
import org.apache.commons.math.ode.sampling.StepHandler;
import org.apache.commons.math.util.FastMath;
@ -189,38 +189,30 @@ public abstract class EmbeddedRungeKuttaIntegrator
/** {@inheritDoc} */
@Override
public double integrate(final ExpandableFirstOrderDifferentialEquations equations,
final double t0, final double[] z0,
final double t, final double[] z)
public void integrate(final ExpandableStatefulODE equations, final double t)
throws MathIllegalStateException, MathIllegalArgumentException {
sanityChecks(equations, t0, z0, t, z);
sanityChecks(equations, t);
setEquations(equations);
resetEvaluations();
final boolean forward = t > t0;
final boolean forward = t > equations.getTime();
// create some internal working arrays
final int totalDim = equations.getDimension();
final int mainDim = equations.getMainSetDimension();
final double[] y0 = new double[totalDim];
final double[] y = new double[totalDim];
System.arraycopy(z0, 0, y0, 0, mainDim);
System.arraycopy(equations.getCurrentAdditionalStates(), 0, y0, mainDim, totalDim - mainDim);
final double[] y0 = equations.getCompleteState();
final double[] y = y0.clone();
final int stages = c.length + 1;
if (y != y0) {
System.arraycopy(y0, 0, y, 0, totalDim);
}
final double[][] yDotK = new double[stages][totalDim];
final double[] yTmp = new double[totalDim];
final double[] yDotTmp = new double[totalDim];
final double[][] yDotK = new double[stages][y.length];
final double[] yTmp = new double[y.length];
final double[] yDotTmp = new double[y.length];
// set up an interpolator sharing the integrator arrays
final RungeKuttaStepInterpolator interpolator = (RungeKuttaStepInterpolator) prototype.copy();
interpolator.reinitialize(this, yTmp, yDotK, forward);
interpolator.storeTime(t0);
interpolator.reinitialize(this, yTmp, yDotK, forward,
equations.getPrimaryMapper(), equations.getSecondaryMappers());
interpolator.storeTime(equations.getTime());
// set up integration control objects
stepStart = t0;
stepStart = equations.getTime();
double hNew = 0;
boolean firstTime = true;
for (StepHandler handler : stepHandlers) {
@ -331,13 +323,11 @@ public abstract class EmbeddedRungeKuttaIntegrator
} while (!isLastStep);
// dispatch result between main and additional states
System.arraycopy(y, 0, z, 0, z.length);
equations.setCurrentAdditionalState(y);
// dispatch results
equations.setTime(stepStart);
equations.setCompleteState(y);
final double stopTime = stepStart;
resetInternalState();
return stopTime;
}

45
src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java
View File

@ -20,7 +20,7 @@ package org.apache.commons.math.ode.nonstiff;
import org.apache.commons.math.analysis.solvers.UnivariateRealSolver;
import org.apache.commons.math.exception.MathIllegalArgumentException;
import org.apache.commons.math.exception.MathIllegalStateException;
import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations;
import org.apache.commons.math.ode.ExpandableStatefulODE;
import org.apache.commons.math.ode.events.EventHandler;
import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
import org.apache.commons.math.ode.sampling.StepHandler;
@ -541,32 +541,27 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
/** {@inheritDoc} */
@Override
public double integrate(final ExpandableFirstOrderDifferentialEquations equations,
final double t0, final double[] z0, final double t, final double[] z)
public void integrate(final ExpandableStatefulODE equations, final double t)
throws MathIllegalStateException, MathIllegalArgumentException {
sanityChecks(equations, t0, z0, t, z);
sanityChecks(equations, t);
setEquations(equations);
resetEvaluations();
final boolean forward = t > t0;
final boolean forward = t > equations.getTime();
// create some internal working arrays
final int totalDim = equations.getDimension();
final int mainDim = equations.getMainSetDimension();
final double[] y0 = new double[totalDim];
final double[] y = new double[totalDim];
System.arraycopy(z0, 0, y0, 0, mainDim);
System.arraycopy(equations.getCurrentAdditionalStates(), 0, y0, mainDim, totalDim - mainDim);
final double[] yDot0 = new double[totalDim];
final double[] y1 = new double[totalDim];
final double[] yTmp = new double[totalDim];
final double[] yTmpDot = new double[totalDim];
final double[] y0 = equations.getCompleteState();
final double[] y = y0.clone();
final double[] yDot0 = new double[y.length];
final double[] y1 = new double[y.length];
final double[] yTmp = new double[y.length];
final double[] yTmpDot = new double[y.length];
final double[][] diagonal = new double[sequence.length-1][];
final double[][] y1Diag = new double[sequence.length-1][];
for (int k = 0; k < sequence.length-1; ++k) {
diagonal[k] = new double[totalDim];
y1Diag[k] = new double[totalDim];
diagonal[k] = new double[y.length];
y1Diag[k] = new double[y.length];
}
final double[][][] fk = new double[sequence.length][][];
@ -606,10 +601,12 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
final AbstractStepInterpolator interpolator =
new GraggBulirschStoerStepInterpolator(y, yDot0,
y1, yDot1,
yMidDots, forward);
interpolator.storeTime(t0);
yMidDots, forward,
equations.getPrimaryMapper(),
equations.getSecondaryMappers());
interpolator.storeTime(equations.getTime());
stepStart = t0;
stepStart = equations.getTime();
double hNew = 0;
double maxError = Double.MAX_VALUE;
boolean previousRejected = false;
@ -940,13 +937,11 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
} while (!isLastStep);
// dispatch result between main and additional states
System.arraycopy(y, 0, z, 0, z.length);
equations.setCurrentAdditionalState(y);
// dispatch results
equations.setTime(stepStart);
equations.setCompleteState(y);
final double stopTime = stepStart;
resetInternalState();
return stopTime;
}

85
src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java
View File

@ -21,6 +21,7 @@ import java.io.IOException;
import java.io.ObjectInput;
import java.io.ObjectOutput;
import org.apache.commons.math.ode.EquationsMapper;
import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
import org.apache.commons.math.ode.sampling.StepInterpolator;
import org.apache.commons.math.util.FastMath;
@ -94,13 +95,13 @@ class GraggBulirschStoerStepInterpolator
*/
private double[][] yMidDots;
/** Interpolation polynoms. */
private double[][] polynoms;
/** Interpolation polynomials. */
private double[][] polynomials;
/** Error coefficients for the interpolation. */
private double[] errfac;
/** Degree of the interpolation polynoms. */
/** Degree of the interpolation polynomials. */
private int currentDegree;
/** Simple constructor.
@ -126,13 +127,17 @@ class GraggBulirschStoerStepInterpolator
* @param yMidDots reference to the integrator array holding the
* derivatives at the middle point of the step
* @param forward integration direction indicator
* @param primaryMapper equations mapper for the primary equations set
* @param secondaryMappers equations mappers for the secondary equations sets
*/
public GraggBulirschStoerStepInterpolator(final double[] y, final double[] y0Dot,
final double[] y1, final double[] y1Dot,
final double[][] yMidDots,
final boolean forward) {
final boolean forward,
final EquationsMapper primaryMapper,
final EquationsMapper[] secondaryMappers) {
super(y, forward);
super(y, forward, primaryMapper, secondaryMappers);
this.y0Dot = y0Dot;
this.y1 = y1;
this.y1Dot = y1Dot;
@ -161,16 +166,16 @@ class GraggBulirschStoerStepInterpolator
y1Dot = null;
yMidDots = null;
// copy the interpolation polynoms (up to the current degree only)
if (interpolator.polynoms == null) {
polynoms = null;
// copy the interpolation polynomials (up to the current degree only)
if (interpolator.polynomials == null) {
polynomials = null;
currentDegree = -1;
} else {
resetTables(interpolator.currentDegree);
for (int i = 0; i < polynoms.length; ++i) {
polynoms[i] = new double[dimension];
System.arraycopy(interpolator.polynoms[i], 0,
polynoms[i], 0, dimension);
for (int i = 0; i < polynomials.length; ++i) {
polynomials[i] = new double[dimension];
System.arraycopy(interpolator.polynomials[i], 0,
polynomials[i], 0, dimension);
}
currentDegree = interpolator.currentDegree;
}
@ -179,21 +184,21 @@ class GraggBulirschStoerStepInterpolator
/** Reallocate the internal tables.
* Reallocate the internal tables in order to be able to handle
* interpolation polynoms up to the given degree
* interpolation polynomials up to the given degree
* @param maxDegree maximal degree to handle
*/
private void resetTables(final int maxDegree) {
if (maxDegree < 0) {
polynoms = null;
polynomials = null;
errfac = null;
currentDegree = -1;
} else {
final double[][] newPols = new double[maxDegree + 1][];
if (polynoms != null) {
System.arraycopy(polynoms, 0, newPols, 0, polynoms.length);
for (int i = polynoms.length; i < newPols.length; ++i) {
if (polynomials != null) {
System.arraycopy(polynomials, 0, newPols, 0, polynomials.length);
for (int i = polynomials.length; i < newPols.length; ++i) {
newPols[i] = new double[currentState.length];
}
} else {
@ -201,7 +206,7 @@ class GraggBulirschStoerStepInterpolator
newPols[i] = new double[currentState.length];
}
}
polynoms = newPols;
polynomials = newPols;
// initialize the error factors array for interpolation
if (maxDegree <= 4) {
@ -237,7 +242,7 @@ class GraggBulirschStoerStepInterpolator
*/
public void computeCoefficients(final int mu, final double h) {
if ((polynoms == null) || (polynoms.length <= (mu + 4))) {
if ((polynomials == null) || (polynomials.length <= (mu + 4))) {
resetTables(mu + 4);
}
@ -251,10 +256,10 @@ class GraggBulirschStoerStepInterpolator
final double aspl = ydiff - yp1;
final double bspl = yp0 - ydiff;
polynoms[0][i] = currentState[i];
polynoms[1][i] = ydiff;
polynoms[2][i] = aspl;
polynoms[3][i] = bspl;
polynomials[0][i] = currentState[i];
polynomials[1][i] = ydiff;
polynomials[2][i] = aspl;
polynomials[3][i] = bspl;
if (mu < 0) {
return;
@ -262,25 +267,25 @@ class GraggBulirschStoerStepInterpolator
// compute the remaining coefficients
final double ph0 = 0.5 * (currentState[i] + y1[i]) + 0.125 * (aspl + bspl);
polynoms[4][i] = 16 * (yMidDots[0][i] - ph0);
polynomials[4][i] = 16 * (yMidDots[0][i] - ph0);
if (mu > 0) {
final double ph1 = ydiff + 0.25 * (aspl - bspl);
polynoms[5][i] = 16 * (yMidDots[1][i] - ph1);
polynomials[5][i] = 16 * (yMidDots[1][i] - ph1);
if (mu > 1) {
final double ph2 = yp1 - yp0;
polynoms[6][i] = 16 * (yMidDots[2][i] - ph2 + polynoms[4][i]);
polynomials[6][i] = 16 * (yMidDots[2][i] - ph2 + polynomials[4][i]);
if (mu > 2) {
final double ph3 = 6 * (bspl - aspl);
polynoms[7][i] = 16 * (yMidDots[3][i] - ph3 + 3 * polynoms[5][i]);
polynomials[7][i] = 16 * (yMidDots[3][i] - ph3 + 3 * polynomials[5][i]);
for (int j = 4; j <= mu; ++j) {
final double fac1 = 0.5 * j * (j - 1);
final double fac2 = 2 * fac1 * (j - 2) * (j - 3);
polynoms[j+4][i] =
16 * (yMidDots[j][i] + fac1 * polynoms[j+2][i] - fac2 * polynoms[j][i]);
polynomials[j+4][i] =
16 * (yMidDots[j][i] + fac1 * polynomials[j+2][i] - fac2 * polynomials[j][i]);
}
}
@ -298,7 +303,7 @@ class GraggBulirschStoerStepInterpolator
double error = 0;
if (currentDegree >= 5) {
for (int i = 0; i < scale.length; ++i) {
final double e = polynoms[currentDegree][i] / scale[i];
final double e = polynomials[currentDegree][i] / scale[i];
error += e * e;
}
error = FastMath.sqrt(error / scale.length) * errfac[currentDegree - 5];
@ -309,7 +314,7 @@ class GraggBulirschStoerStepInterpolator
/** {@inheritDoc} */
@Override
protected void computeInterpolatedStateAndDerivatives(final double theta,
final double oneMinusThetaH) {
final double oneMinusThetaH) {
final int dimension = currentState.length;
@ -324,20 +329,20 @@ class GraggBulirschStoerStepInterpolator
for (int i = 0; i < dimension; ++i) {
final double p0 = polynoms[0][i];
final double p1 = polynoms[1][i];
final double p2 = polynoms[2][i];
final double p3 = polynoms[3][i];
final double p0 = polynomials[0][i];
final double p1 = polynomials[1][i];
final double p2 = polynomials[2][i];
final double p3 = polynomials[3][i];
interpolatedState[i] = p0 + theta * (p1 + oneMinusTheta * (p2 * theta + p3 * oneMinusTheta));
interpolatedDerivatives[i] = dot1 * p1 + dot2 * p2 + dot3 * p3;
if (currentDegree > 3) {
double cDot = 0;
double c = polynoms[currentDegree][i];
double c = polynomials[currentDegree][i];
for (int j = currentDegree - 1; j > 3; --j) {
final double d = 1.0 / (j - 3);
cDot = d * (theta05 * cDot + c);
c = polynoms[j][i] + c * d * theta05;
c = polynomials[j][i] + c * d * theta05;
}
interpolatedState[i] += t4 * c;
interpolatedDerivatives[i] += (t4 * cDot + t4Dot * c) / h;
@ -367,7 +372,7 @@ class GraggBulirschStoerStepInterpolator
out.writeInt(currentDegree);
for (int k = 0; k <= currentDegree; ++k) {
for (int l = 0; l < dimension; ++l) {
out.writeDouble(polynoms[k][l]);
out.writeDouble(polynomials[k][l]);
}
}
@ -376,7 +381,7 @@ class GraggBulirschStoerStepInterpolator
/** {@inheritDoc} */
@Override
public void readExternal(final ObjectInput in)
throws IOException {
throws IOException, ClassNotFoundException {
// read the base class
final double t = readBaseExternal(in);
@ -389,7 +394,7 @@ class GraggBulirschStoerStepInterpolator
for (int k = 0; k <= currentDegree; ++k) {
for (int l = 0; l < dimension; ++l) {
polynoms[k][l] = in.readDouble();
polynomials[k][l] = in.readDouble();
}
}

38
src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java
View File

@ -21,7 +21,7 @@ package org.apache.commons.math.ode.nonstiff;
import org.apache.commons.math.exception.MathIllegalArgumentException;
import org.apache.commons.math.exception.MathIllegalStateException;
import org.apache.commons.math.ode.AbstractIntegrator;
import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations;
import org.apache.commons.math.ode.ExpandableStatefulODE;
import org.apache.commons.math.ode.sampling.StepHandler;
import org.apache.commons.math.util.FastMath;
@ -90,27 +90,18 @@ public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
}
/** {@inheritDoc} */
public double integrate(final ExpandableFirstOrderDifferentialEquations equations,
final double t0, final double[] z0,
final double t, final double[] z)
public void integrate(final ExpandableStatefulODE equations, final double t)
throws MathIllegalStateException, MathIllegalArgumentException {
sanityChecks(equations, t0, z0, t, z);
sanityChecks(equations, t);
setEquations(equations);
resetEvaluations();
final boolean forward = t > t0;
final boolean forward = t > equations.getTime();
// create some internal working arrays
final int totalDim = equations.getDimension();
final int mainDim = equations.getMainSetDimension();
final double[] y0 = new double[totalDim];
final double[] y = new double[totalDim];
System.arraycopy(z0, 0, y0, 0, mainDim);
System.arraycopy(equations.getCurrentAdditionalStates(), 0, y0, mainDim, totalDim - mainDim);
final int stages = c.length + 1;
if (y != y0) {
System.arraycopy(y0, 0, y, 0, y0.length);
}
final double[] y0 = equations.getCompleteState();
final double[] y = y0.clone();
final int stages = c.length + 1;
final double[][] yDotK = new double[stages][];
for (int i = 0; i < stages; ++i) {
yDotK [i] = new double[y0.length];
@ -120,11 +111,12 @@ public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
// set up an interpolator sharing the integrator arrays
final RungeKuttaStepInterpolator interpolator = (RungeKuttaStepInterpolator) prototype.copy();
interpolator.reinitialize(this, yTmp, yDotK, forward);
interpolator.storeTime(t0);
interpolator.reinitialize(this, yTmp, yDotK, forward,
equations.getPrimaryMapper(), equations.getSecondaryMappers());
interpolator.storeTime(equations.getTime());
// set up integration control objects
stepStart = t0;
stepStart = equations.getTime();
stepSize = forward ? step : -step;
for (StepHandler handler : stepHandlers) {
handler.reset();
@ -185,14 +177,12 @@ public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
} while (!isLastStep);
// dispatch result between main and additional states
System.arraycopy(y, 0, z, 0, z.length);
equations.setCurrentAdditionalState(y);
// dispatch results
equations.setTime(stepStart);
equations.setCompleteState(y);
final double stopTime = stepStart;
stepStart = Double.NaN;
stepSize = Double.NaN;
return stopTime;
}

11
src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaStepInterpolator.java
View File

@ -22,6 +22,7 @@ import java.io.ObjectInput;
import java.io.ObjectOutput;
import org.apache.commons.math.ode.AbstractIntegrator;
import org.apache.commons.math.ode.EquationsMapper;
import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
/** This class represents an interpolator over the last step during an
@ -120,10 +121,14 @@ abstract class RungeKuttaStepInterpolator
* @param yDotArray reference to the integrator array holding all the
* intermediate slopes
* @param forward integration direction indicator
* @param primaryMapper equations mapper for the primary equations set
* @param secondaryMappers equations mappers for the secondary equations sets
*/
public void reinitialize(final AbstractIntegrator rkIntegrator,
final double[] y, final double[][] yDotArray, final boolean forward) {
reinitialize(y, forward);
final double[] y, final double[][] yDotArray, final boolean forward,
final EquationsMapper primaryMapper,
final EquationsMapper[] secondaryMappers) {
reinitialize(y, forward, primaryMapper, secondaryMappers);
this.yDotK = yDotArray;
this.integrator = rkIntegrator;
}
@ -153,7 +158,7 @@ abstract class RungeKuttaStepInterpolator
/** {@inheritDoc} */
@Override
public void readExternal(final ObjectInput in)
throws IOException {
throws IOException, ClassNotFoundException {
// read the base class
final double t = readBaseExternal(in);

240
src/main/java/org/apache/commons/math/ode/sampling/AbstractStepInterpolator.java
View File

@ -21,6 +21,8 @@ import java.io.IOException;
import java.io.ObjectInput;
import java.io.ObjectOutput;
import org.apache.commons.math.ode.EquationsMapper;
/** This abstract class represents an interpolator over the last step
* during an ODE integration.
*
@ -56,6 +58,18 @@ public abstract class AbstractStepInterpolator
/** interpolated derivatives */
protected double[] interpolatedDerivatives;
/** interpolated primary state */
protected double[] interpolatedPrimaryState;
/** interpolated primary derivatives */
protected double[] interpolatedPrimaryDerivatives;
/** interpolated secondary state */
protected double[][] interpolatedSecondaryState;
/** interpolated secondary derivatives */
protected double[][] interpolatedSecondaryDerivatives;
/** global previous time */
private double globalPreviousTime;
@ -77,6 +91,11 @@ public abstract class AbstractStepInterpolator
/** indicator for dirty state. */
private boolean dirtyState;
/** Equations mapper for the primary equations set. */
private EquationsMapper primaryMapper;
/** Equations mappers for the secondary equations sets. */
private EquationsMapper[] secondaryMappers;
/** Simple constructor.
* This constructor builds an instance that is not usable yet, the
@ -90,41 +109,45 @@ public abstract class AbstractStepInterpolator
* initializing the copy.
*/
protected AbstractStepInterpolator() {
globalPreviousTime = Double.NaN;
globalCurrentTime = Double.NaN;
softPreviousTime = Double.NaN;
softCurrentTime = Double.NaN;
h = Double.NaN;
interpolatedTime = Double.NaN;
currentState = null;
interpolatedState = null;
interpolatedDerivatives = null;
finalized = false;
this.forward = true;
this.dirtyState = true;
}
/** Simple constructor.
* @param y reference to the integrator array holding the state at
* the end of the step
* @param forward integration direction indicator
*/
protected AbstractStepInterpolator(final double[] y, final boolean forward) {
globalPreviousTime = Double.NaN;
globalCurrentTime = Double.NaN;
softPreviousTime = Double.NaN;
softCurrentTime = Double.NaN;
h = Double.NaN;
interpolatedTime = Double.NaN;
currentState = null;
finalized = false;
this.forward = true;
this.dirtyState = true;
primaryMapper = null;
secondaryMappers = null;
allocateInterpolatedArrays(-1, null, null);
}
currentState = y;
interpolatedState = new double[y.length];
interpolatedDerivatives = new double[y.length];
/** Simple constructor.
* @param y reference to the integrator array holding the state at
* the end of the step
* @param forward integration direction indicator
* @param primaryMapper equations mapper for the primary equations set
* @param secondaryMappers equations mappers for the secondary equations sets
*/
protected AbstractStepInterpolator(final double[] y, final boolean forward,
final EquationsMapper primaryMapper,
final EquationsMapper[] secondaryMappers) {
finalized = false;
this.forward = forward;
this.dirtyState = true;
globalPreviousTime = Double.NaN;
globalCurrentTime = Double.NaN;
softPreviousTime = Double.NaN;
softCurrentTime = Double.NaN;
h = Double.NaN;
interpolatedTime = Double.NaN;
currentState = y;
finalized = false;
this.forward = forward;
this.dirtyState = true;
this.primaryMapper = primaryMapper;
this.secondaryMappers = (secondaryMappers == null) ? null : secondaryMappers.clone();
allocateInterpolatedArrays(y.length, primaryMapper, secondaryMappers);
}
@ -154,43 +177,89 @@ public abstract class AbstractStepInterpolator
h = interpolator.h;
interpolatedTime = interpolator.interpolatedTime;
if (interpolator.currentState != null) {
currentState = interpolator.currentState.clone();
interpolatedState = interpolator.interpolatedState.clone();
interpolatedDerivatives = interpolator.interpolatedDerivatives.clone();
if (interpolator.currentState == null) {
currentState = null;
allocateInterpolatedArrays(-1, null, null);
} else {
currentState = null;
interpolatedState = null;
interpolatedDerivatives = null;
currentState = interpolator.currentState.clone();
interpolatedState = interpolator.interpolatedState.clone();
interpolatedDerivatives = interpolator.interpolatedDerivatives.clone();
interpolatedPrimaryState = interpolator.interpolatedPrimaryState.clone();
interpolatedPrimaryDerivatives = interpolator.interpolatedPrimaryDerivatives.clone();
interpolatedSecondaryState = new double[interpolator.interpolatedSecondaryState.length][];
interpolatedSecondaryDerivatives = new double[interpolator.interpolatedSecondaryDerivatives.length][];
for (int i = 0; i < interpolatedSecondaryState.length; ++i) {
interpolatedSecondaryState[i] = interpolator.interpolatedSecondaryState[i].clone();
interpolatedSecondaryDerivatives[i] = interpolator.interpolatedSecondaryDerivatives[i].clone();
}
}
finalized = interpolator.finalized;
forward = interpolator.forward;
dirtyState = interpolator.dirtyState;
finalized = interpolator.finalized;
forward = interpolator.forward;
dirtyState = interpolator.dirtyState;
primaryMapper = interpolator.primaryMapper;
secondaryMappers = (interpolator.secondaryMappers == null) ?
null : interpolator.secondaryMappers.clone();
}
/** Reinitialize the instance
* @param y reference to the integrator array holding the state at
* the end of the step
* @param isForward integration direction indicator
/** Allocate the various interpolated states arrays.
* @param dimension total dimension (negative if arrays should be set to null)
* @param primaryMapper equations mapper for the primary equations set
* @param secondaryMappers equations mappers for the secondary equations sets
*/
protected void reinitialize(final double[] y, final boolean isForward) {
private void allocateInterpolatedArrays(final int dimension,
final EquationsMapper primaryMapper,
final EquationsMapper[] secondaryMappers) {
if (dimension < 0) {
interpolatedState = null;
interpolatedDerivatives = null;
interpolatedPrimaryState = null;
interpolatedPrimaryDerivatives = null;
interpolatedSecondaryState = null;
interpolatedSecondaryDerivatives = null;
} else {
interpolatedState = new double[dimension];
interpolatedDerivatives = new double[dimension];
interpolatedPrimaryState = new double[primaryMapper.getDimension()];
interpolatedPrimaryDerivatives = new double[primaryMapper.getDimension()];
if (secondaryMappers == null) {
interpolatedSecondaryState = null;
interpolatedSecondaryDerivatives = null;
} else {
interpolatedSecondaryState = new double[secondaryMappers.length][];
interpolatedSecondaryDerivatives = new double[secondaryMappers.length][];
for (int i = 0; i < secondaryMappers.length; ++i) {
interpolatedSecondaryState[i] = new double[secondaryMappers[i].getDimension()];
interpolatedSecondaryDerivatives[i] = new double[secondaryMappers[i].getDimension()];
}
}
}
}
globalPreviousTime = Double.NaN;
globalCurrentTime = Double.NaN;
softPreviousTime = Double.NaN;
softCurrentTime = Double.NaN;
h = Double.NaN;
interpolatedTime = Double.NaN;
/** Reinitialize the instance
* @param y reference to the integrator array holding the state at the end of the step
* @param isForward integration direction indicator
* @param primaryMapper equations mapper for the primary equations set
* @param secondaryMappers equations mappers for the secondary equations sets
*/
protected void reinitialize(final double[] y, final boolean isForward,
final EquationsMapper primaryMapper,
final EquationsMapper[] secondaryMappers) {
currentState = y;
interpolatedState = new double[y.length];
interpolatedDerivatives = new double[y.length];
finalized = false;
this.forward = isForward;
this.dirtyState = true;
globalPreviousTime = Double.NaN;
globalCurrentTime = Double.NaN;
softPreviousTime = Double.NaN;
softCurrentTime = Double.NaN;
h = Double.NaN;
interpolatedTime = Double.NaN;
currentState = y;
finalized = false;
this.forward = isForward;
this.dirtyState = true;
this.primaryMapper = primaryMapper;
this.secondaryMappers = secondaryMappers.clone();
allocateInterpolatedArrays(y.length, primaryMapper, secondaryMappers);
}
@ -328,9 +397,9 @@ public abstract class AbstractStepInterpolator
protected abstract void computeInterpolatedStateAndDerivatives(double theta,
double oneMinusThetaH);
/** {@inheritDoc} */
public double[] getInterpolatedState() {
/** Lazy evaluation of complete interpolated state.
*/
private void evaluateCompleteInterpolatedState() {
// lazy evaluation of the state
if (dirtyState) {
final double oneMinusThetaH = globalCurrentTime - interpolatedTime;
@ -338,24 +407,38 @@ public abstract class AbstractStepInterpolator
computeInterpolatedStateAndDerivatives(theta, oneMinusThetaH);
dirtyState = false;
}
}
return interpolatedState;
/** {@inheritDoc} */
public double[] getInterpolatedState() {
evaluateCompleteInterpolatedState();
primaryMapper.extractEquationData(interpolatedState,
interpolatedPrimaryState);
return interpolatedPrimaryState;
}
/** {@inheritDoc} */
public double[] getInterpolatedDerivatives() {
evaluateCompleteInterpolatedState();
primaryMapper.extractEquationData(interpolatedDerivatives,
interpolatedPrimaryDerivatives);
return interpolatedPrimaryDerivatives;
}
// lazy evaluation of the state
if (dirtyState) {
final double oneMinusThetaH = globalCurrentTime - interpolatedTime;
final double theta = (h == 0) ? 0 : (h - oneMinusThetaH) / h;
computeInterpolatedStateAndDerivatives(theta, oneMinusThetaH);
dirtyState = false;
}
return interpolatedDerivatives;
/** {@inheritDoc} */
public double[] getInterpolatedSecondaryState(final int index) {
evaluateCompleteInterpolatedState();
secondaryMappers[index].extractEquationData(interpolatedState,
interpolatedSecondaryState[index]);
return interpolatedSecondaryState[index];
}
/** {@inheritDoc} */
public double[] getInterpolatedSecondaryDerivatives(final int index) {
evaluateCompleteInterpolatedState();
secondaryMappers[index].extractEquationData(interpolatedDerivatives,
interpolatedSecondaryDerivatives[index]);
return interpolatedSecondaryDerivatives[index];
}
/**
@ -439,6 +522,11 @@ public abstract class AbstractStepInterpolator
out.writeDouble(softCurrentTime);
out.writeDouble(h);
out.writeBoolean(forward);
out.writeObject(primaryMapper);
out.write(secondaryMappers.length);
for (final EquationsMapper mapper : secondaryMappers) {
out.writeObject(mapper);
}
if (currentState != null) {
for (int i = 0; i < currentState.length; ++i) {
@ -468,11 +556,11 @@ public abstract class AbstractStepInterpolator
* properly calling the {@link #setInterpolatedTime} method later,
* once all rest of the object state has been set up properly.
* @param in stream where to read the state from
* @return interpolated time be set later by the caller
* @return interpolated time to be set later by the caller
* @exception IOException in case of read error
*/
protected double readBaseExternal(final ObjectInput in)
throws IOException {
throws IOException, ClassNotFoundException {
final int dimension = in.readInt();
globalPreviousTime = in.readDouble();
@ -481,6 +569,11 @@ public abstract class AbstractStepInterpolator
softCurrentTime = in.readDouble();
h = in.readDouble();
forward = in.readBoolean();
primaryMapper = (EquationsMapper) in.readObject();
secondaryMappers = new EquationsMapper[in.read()];
for (int i = 0; i < secondaryMappers.length; ++i) {
secondaryMappers[i] = (EquationsMapper) in.readObject();
}
dirtyState = true;
if (dimension < 0) {
@ -493,9 +586,8 @@ public abstract class AbstractStepInterpolator
}
// we do NOT handle the interpolated time and state here
interpolatedTime = Double.NaN;
interpolatedState = (dimension < 0) ? null : new double[dimension];
interpolatedDerivatives = (dimension < 0) ? null : new double[dimension];
interpolatedTime = Double.NaN;
allocateInterpolatedArrays(dimension, primaryMapper, secondaryMappers);
finalized = true;

9
src/main/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolator.java
View File

@ -23,6 +23,7 @@ import java.io.ObjectOutput;
import java.util.Arrays;
import org.apache.commons.math.linear.Array2DRowRealMatrix;
import org.apache.commons.math.ode.EquationsMapper;
import org.apache.commons.math.util.FastMath;
/**
@ -104,10 +105,14 @@ public class NordsieckStepInterpolator extends AbstractStepInterpolator {
* @param y reference to the integrator array holding the state at
* the end of the step
* @param forward integration direction indicator
* @param primaryMapper equations mapper for the primary equations set
* @param secondaryMappers equations mappers for the secondary equations sets
*/
@Override
public void reinitialize(final double[] y, final boolean forward) {
super.reinitialize(y, forward);
public void reinitialize(final double[] y, final boolean forward,
final EquationsMapper primaryMapper,
final EquationsMapper[] secondaryMappers) {
super.reinitialize(y, forward, primaryMapper, secondaryMappers);
stateVariation = new double[y.length];
}

34
src/main/java/org/apache/commons/math/ode/sampling/StepInterpolator.java
View File

@ -101,6 +101,40 @@ public interface StepInterpolator extends Externalizable {
*/
double[] getInterpolatedDerivatives();
/** Get the interpolated secondary state corresponding to the secondary equations.
* <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>
* @param index index of the secondary set, as returned by {@link
* org.apache.commons.math.ode.ExpandableStatefulODE#addSecondaryEquations(
* org.apache.commons.math.ode.SecondaryEquations)
* ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)}
* @return interpolated secondary state at the current interpolation date
* @see #getInterpolatedState()
* @see #getInterpolatedDerivatives()
* @see #getInterpolatedSecondaryDerivatives(String)
* @see #setInterpolatedTime(double)
* @since 3.0
*/
double[] getInterpolatedSecondaryState(int index);
/** Get the interpolated secondary derivatives corresponding to the secondary equations.
* <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>
* @param index index of the secondary set, as returned by {@link
* org.apache.commons.math.ode.ExpandableStatefulODE#addSecondaryEquations(
* org.apache.commons.math.ode.SecondaryEquations)
* ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)}
* @return interpolated secondary derivatives at the current interpolation date
* @see #getInterpolatedState()
* @see #getInterpolatedDerivatives()
* @see #getInterpolatedSecondaryState(String)
* @see #setInterpolatedTime(double)
* @since 3.0
*/
double[] getInterpolatedSecondaryDerivatives(int index);
/** 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

2
src/main/resources/META-INF/localization/LocalizedFormats_fr.properties
View File

@ -281,7 +281,7 @@ UNBOUNDED_SOLUTION = solution non born\u00e9e
UNKNOWN_MODE = mode {0} inconnu, modes connus : {1} ({2}), {3} ({4}), {5} ({6}), {7} ({8}), {9} ({10}) et {11} ({12})
UNKNOWN_ADDITIONAL_EQUATION = \u00e9quation additionnelle inconnue
UNKNOWN_PARAMETER = param\u00e8tre {0} inconnu
UNMATCHED_ODE_IN_EXTENDED_SET = l''\u00e9quation diff\u00e9rentielle ne correspond pas \u00e0 l''\u00e9quation principale du jeu \u00e9tendu
UNMATCHED_ODE_IN_EXPANDED_SET = l''\u00e9quation diff\u00e9rentielle ne correspond pas \u00e0 l''\u00e9quation principale du jeu \u00e9tendu
CANNOT_PARSE_AS_TYPE = cha\u00eene {0} non analysable (\u00e0 partir de la position {1}) en un objet de type {2}
CANNOT_PARSE = cha\u00eene {0} non analysable (\u00e0 partir de la position {1})
UNPARSEABLE_3D_VECTOR = vecteur 3D non analysable : "{0}"

70
src/test/java/org/apache/commons/math/ode/JacobianMatricesTest.java
View File

@ -85,7 +85,7 @@ public class JacobianMatricesTest {
@Test
public void testInternalDifferentiation() {
ExpandableFirstOrderIntegrator integ =
AbstractIntegrator 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;
double hY = 1.0e-12;
@ -97,15 +97,19 @@ public class JacobianMatricesTest {
double[] z = { 1.3, b };
double[][] dZdZ0 = new double[2][2];
double[] dZdP = new double[2];
ExpandableFirstOrderDifferentialEquations efode = new ExpandableFirstOrderDifferentialEquations(brusselator);
JacobianMatrices jacob = new JacobianMatrices(efode, brusselator);
JacobianMatrices jacob = new JacobianMatrices(brusselator, new double[] { hY, hY }, ParamBrusselator.B);
jacob.setParameterizedODE(brusselator);
jacob.selectParameters(ParamBrusselator.B);
jacob.setMainStateSteps(new double[] { hY, hY });
jacob.setParameterStep(ParamBrusselator.B, hP);
jacob.setInitialParameterJacobian(ParamBrusselator.B, new double[] { 0.0, 1.0 });
ExpandableStatefulODE efode = new ExpandableStatefulODE(brusselator);
efode.setTime(0);
efode.setPrimaryState(z);
jacob.registerVariationalEquations(efode);
integ.setMaxEvaluations(5000);
integ.integrate(efode, 0, z, 20.0, z);
integ.integrate(efode, 20.0);
jacob.getCurrentMainSetJacobian(dZdZ0);
jacob.getCurrentParameterJacobian(ParamBrusselator.B, dZdP);
// Assert.assertEquals(5000, integ.getMaxEvaluations());
@ -123,7 +127,7 @@ public class JacobianMatricesTest {
@Test
public void testAnalyticalDifferentiation() {
ExpandableFirstOrderIntegrator integ =
AbstractIntegrator 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();
@ -132,13 +136,18 @@ public class JacobianMatricesTest {
double[] z = { 1.3, b };
double[][] dZdZ0 = new double[2][2];
double[] dZdP = new double[2];
ExpandableFirstOrderDifferentialEquations efode = new ExpandableFirstOrderDifferentialEquations(brusselator);
JacobianMatrices jacob = new JacobianMatrices(efode, brusselator);
jacob.setParameterJacobianProvider(brusselator);
jacob.selectParameters(Brusselator.B);
JacobianMatrices jacob = new JacobianMatrices(brusselator, Brusselator.B);
jacob.addParameterJacobianProvider(brusselator);
jacob.setInitialParameterJacobian(Brusselator.B, new double[] { 0.0, 1.0 });
ExpandableStatefulODE efode = new ExpandableStatefulODE(brusselator);
efode.setTime(0);
efode.setPrimaryState(z);
jacob.registerVariationalEquations(efode);
integ.setMaxEvaluations(5000);
integ.integrate(efode, 0, z, 20.0, z);
integ.integrate(efode, 20.0);
jacob.getCurrentMainSetJacobian(dZdZ0);
jacob.getCurrentParameterJacobian(Brusselator.B, dZdP);
// Assert.assertEquals(5000, integ.getMaxEvaluations());
@ -156,23 +165,28 @@ public class JacobianMatricesTest {
@Test
public void testFinalResult() {
ExpandableFirstOrderIntegrator integ =
AbstractIntegrator 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);
ExpandableFirstOrderDifferentialEquations efode = new ExpandableFirstOrderDifferentialEquations(circle);
JacobianMatrices jacob = new JacobianMatrices(efode, circle);
jacob.setParameterJacobianProvider(circle);
jacob.selectParameters(Circle.CX, Circle.CY, Circle.OMEGA);
JacobianMatrices jacob = new JacobianMatrices(circle, Circle.CX, Circle.CY, Circle.OMEGA);
jacob.addParameterJacobianProvider(circle);
jacob.setInitialMainStateJacobian(circle.exactDyDy0(0));
jacob.setInitialParameterJacobian(Circle.CX, circle.exactDyDcx(0));
jacob.setInitialParameterJacobian(Circle.CY, circle.exactDyDcy(0));
jacob.setInitialParameterJacobian(Circle.OMEGA, circle.exactDyDom(0));
ExpandableStatefulODE efode = new ExpandableStatefulODE(circle);
efode.setTime(0);
efode.setPrimaryState(y);
jacob.registerVariationalEquations(efode);
integ.setMaxEvaluations(5000);
double t = 18 * FastMath.PI;
integ.integrate(efode, 0, y, t, y);
integ.integrate(efode, t);
y = efode.getPrimaryState();
for (int i = 0; i < y.length; ++i) {
Assert.assertEquals(circle.exactY(t)[i], y[i], 1.0e-9);
}
@ -204,7 +218,7 @@ public class JacobianMatricesTest {
@Test
public void testParameterizable() {
ExpandableFirstOrderIntegrator integ =
AbstractIntegrator 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 };
ParameterizedCircle pcircle = new ParameterizedCircle(y, 1.0, 1.0, 0.1);
@ -212,25 +226,29 @@ public class JacobianMatricesTest {
// pcircle.setParameter(ParameterizedCircle.CY, 1.0);
// pcircle.setParameter(ParameterizedCircle.OMEGA, 0.1);
ExpandableFirstOrderDifferentialEquations efode = new ExpandableFirstOrderDifferentialEquations(pcircle);
double hP = 1.0e-12;
double hY = 1.0e-12;
JacobianMatrices jacob = new JacobianMatrices(efode, pcircle);
jacob.setParameterJacobianProvider(pcircle);
JacobianMatrices jacob = new JacobianMatrices(pcircle, new double[] { hY, hY },
Circle.CX, Circle.OMEGA);
jacob.addParameterJacobianProvider(pcircle);
jacob.setParameterizedODE(pcircle);
jacob.selectParameters(Circle.CX, Circle.OMEGA);
jacob.setMainStateSteps(new double[] { hY, hY });
jacob.setParameterStep(Circle.OMEGA, hP);
jacob.setInitialMainStateJacobian(pcircle.exactDyDy0(0));
jacob.setInitialParameterJacobian(Circle.CX, pcircle.exactDyDcx(0));
// jacob.setInitialParameterJacobian(Circle.CY, circle.exactDyDcy(0));
jacob.setInitialParameterJacobian(Circle.OMEGA, pcircle.exactDyDom(0));
ExpandableStatefulODE efode = new ExpandableStatefulODE(pcircle);
efode.setTime(0);
efode.setPrimaryState(y);
jacob.registerVariationalEquations(efode);
integ.setMaxEvaluations(50000);
double t = 18 * FastMath.PI;
integ.integrate(efode, 0, y, t, y);
integ.integrate(efode, t);
y = efode.getPrimaryState();
for (int i = 0; i < y.length; ++i) {
Assert.assertEquals(pcircle.exactY(t)[i], y[i], 1.0e-9);
}

13
src/test/java/org/apache/commons/math/ode/nonstiff/EulerStepInterpolatorTest.java
View File

@ -26,6 +26,7 @@ import java.io.ObjectOutputStream;
import java.util.Random;
import org.apache.commons.math.ode.ContinuousOutputModel;
import org.apache.commons.math.ode.EquationsMapper;
import org.apache.commons.math.ode.TestProblem1;
import org.apache.commons.math.ode.TestProblem3;
import org.apache.commons.math.ode.sampling.StepHandler;
@ -42,7 +43,9 @@ public class EulerStepInterpolatorTest {
double[] y = { 0.0, 1.0, -2.0 };
double[][] yDot = { { 1.0, 2.0, -2.0 } };
EulerStepInterpolator interpolator = new EulerStepInterpolator();
interpolator.reinitialize(new DummyIntegrator(interpolator), y, yDot, true);
interpolator.reinitialize(new DummyIntegrator(interpolator), y, yDot, true,
new EquationsMapper(0, y.length),
new EquationsMapper[0]);
interpolator.storeTime(0);
interpolator.shift();
interpolator.storeTime(1);
@ -63,7 +66,9 @@ public class EulerStepInterpolatorTest {
double[] y = y0.clone();
double[][] yDot = { new double[y0.length] };
EulerStepInterpolator interpolator = new EulerStepInterpolator();
interpolator.reinitialize(new DummyIntegrator(interpolator), y, yDot, true);
interpolator.reinitialize(new DummyIntegrator(interpolator), y, yDot, true,
new EquationsMapper(0, y.length),
new EquationsMapper[0]);
interpolator.storeTime(t0);
double dt = 1.0;
@ -96,7 +101,9 @@ public class EulerStepInterpolatorTest {
double[] y = { 1.0, 3.0, -4.0 };
double[][] yDot = { { 1.0, 2.0, -2.0 } };
EulerStepInterpolator interpolator = new EulerStepInterpolator();
interpolator.reinitialize(new DummyIntegrator(interpolator), y, yDot, true);
interpolator.reinitialize(new DummyIntegrator(interpolator), y, yDot, true,
new EquationsMapper(0, y.length),
new EquationsMapper[0]);
interpolator.storeTime(0);
interpolator.shift();
interpolator.storeTime(1);

8
src/main/java/org/apache/commons/math/ode/sampling/DummyStepInterpolator.java → src/test/java/org/apache/commons/math/ode/sampling/DummyStepInterpolator.java
View File

@ -21,6 +21,8 @@ import java.io.IOException;
import java.io.ObjectInput;
import java.io.ObjectOutput;
import org.apache.commons.math.ode.EquationsMapper;
/** This class is a step interpolator that does nothing.
*
* <p>This class is used when the {@link StepHandler "step handler"}
@ -65,9 +67,11 @@ public class DummyStepInterpolator
* @param yDot reference to the integrator array holding the state
* derivative at some arbitrary point within the step
* @param forward integration direction indicator
* @param primaryMapper equations mapper for the primary equations set
* @param secondaryMappers equations mappers for the secondary equations sets
*/
public DummyStepInterpolator(final double[] y, final double[] yDot, final boolean forward) {
super(y, forward);
super(y, forward, new EquationsMapper(0, y.length), new EquationsMapper[0]);
currentDerivative = yDot;
}
@ -128,7 +132,7 @@ public class DummyStepInterpolator
*/
@Override
public void readExternal(final ObjectInput in)
throws IOException {
throws IOException, ClassNotFoundException {
// read the base class
final double t = readBaseExternal(in);