diff --git a/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java b/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java index 1c766144e..afa3376f4 100644 --- a/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java +++ b/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java @@ -315,6 +315,9 @@ 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"), 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}\""), diff --git a/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java b/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java index ad423c995..ebad65b2e 100644 --- a/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java +++ b/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java @@ -29,6 +29,7 @@ import java.util.TreeSet; import org.apache.commons.math.analysis.solvers.BracketingNthOrderBrentSolver; import org.apache.commons.math.analysis.solvers.UnivariateRealSolver; import org.apache.commons.math.exception.DimensionMismatchException; +import org.apache.commons.math.exception.MathIllegalArgumentException; import org.apache.commons.math.exception.MathIllegalStateException; import org.apache.commons.math.exception.MaxCountExceededException; import org.apache.commons.math.exception.NumberIsTooSmallException; @@ -46,7 +47,7 @@ import org.apache.commons.math.util.MathUtils; * @version $Id$ * @since 2.0 */ -public abstract class AbstractIntegrator implements FirstOrderIntegrator { +public abstract class AbstractIntegrator implements ExpandableFirstOrderIntegrator { /** Step handler. */ protected Collection stepHandlers; @@ -76,7 +77,7 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator { private Incrementor evaluations; /** Differential equations to integrate. */ - private transient FirstOrderDifferentialEquations equations; + private transient ExpandableFirstOrderDifferentialEquations equations; /** Build an instance. * @param name name of the method @@ -188,10 +189,17 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator { * @param equations differential equations to integrate * @see #computeDerivatives(double, double[], double[]) */ - protected void setEquations(final FirstOrderDifferentialEquations equations) { + protected void setEquations(final ExpandableFirstOrderDifferentialEquations equations) { this.equations = equations; } + /** {@inheritDoc} */ + public double integrate(FirstOrderDifferentialEquations equations, + double t0, double[] y0, double t, double[] y) + throws MathIllegalStateException, MathIllegalArgumentException { + return integrate(new ExpandableFirstOrderDifferentialEquations(equations), t0, y0, t, y); + } + /** Compute the derivatives and check the number of evaluations. * @param t current value of the independent time variable * @param y array containing the current value of the state vector @@ -336,16 +344,16 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator { * @exception DimensionMismatchException if some inconsistency is detected * @exception NumberIsTooSmallException if integration span is too small */ - protected void sanityChecks(final FirstOrderDifferentialEquations ode, + protected void sanityChecks(final ExpandableFirstOrderDifferentialEquations ode, final double t0, final double[] y0, final double t, final double[] y) throws DimensionMismatchException, NumberIsTooSmallException { - if (ode.getDimension() != y0.length) { + if (ode.getMainSetDimension() != y0.length) { throw new DimensionMismatchException(ode.getDimension(), y0.length); } - if (ode.getDimension() != y.length) { + if (ode.getMainSetDimension() != y.length) { throw new DimensionMismatchException(ode.getDimension(), y.length); } diff --git a/src/main/java/org/apache/commons/math/ode/AbstractParameterizable.java b/src/main/java/org/apache/commons/math/ode/AbstractParameterizable.java new file mode 100644 index 000000000..4766aa23e --- /dev/null +++ b/src/main/java/org/apache/commons/math/ode/AbstractParameterizable.java @@ -0,0 +1,81 @@ +/* + * 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.Collection; + +import org.apache.commons.math.exception.MathIllegalArgumentException; +import org.apache.commons.math.exception.util.LocalizedFormats; + +/** This abstract class provides boilerplate parameters list. + * + * @version $Id$ + * @since 3.0 + */ + +public abstract class AbstractParameterizable implements Parameterizable { + + /** List of the parameters names. */ + private final Collection parametersNames; + + /** Simple constructor. + * @param names names of the supported parameters + */ + protected AbstractParameterizable(final String ... names) { + parametersNames = new ArrayList(); + for (final String name : names) { + parametersNames.add(name); + } + } + + /** Simple constructor. + * @param names names of the supported parameters + */ + protected AbstractParameterizable(final Collection names) { + parametersNames = new ArrayList(); + parametersNames.addAll(names); + } + + /** {@inheritDoc} */ + public Collection getParametersNames() { + return parametersNames; + } + + /** {@inheritDoc} */ + public boolean isSupported(final String name) { + for (final String supportedName : parametersNames) { + if (supportedName.equals(name)) { + return true; + } + } + return false; + } + + /** Check if a parameter is supported and throw an IllegalArgumentException if not. + * @param name name of the parameter to check + * @exception MathIllegalArgumentException if the parameter is not supported + * @see #isSupported(String) + */ + public void complainIfNotSupported(final String name) + throws MathIllegalArgumentException { + if (!isSupported(name)) { + throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER, name); + } + } + +} diff --git a/src/main/java/org/apache/commons/math/ode/AdditionalEquations.java b/src/main/java/org/apache/commons/math/ode/AdditionalEquations.java new file mode 100644 index 000000000..046aaeb17 --- /dev/null +++ b/src/main/java/org/apache/commons/math/ode/AdditionalEquations.java @@ -0,0 +1,55 @@ +/* + * 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 allows users to add their own differential equations to a main + * set of differential equations. + *

+ * 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 + * adjoined parameters linked to the minimized hamiltonian must be integrated. + *

+ *

+ * This interface allows users to add such equations to a main set of {@link + * FirstOrderDifferentialEquations first order differential equations} + * thanks to the {@link + * ExpandableFirstOrderDifferentialEquations#addAdditionalEquations(AdditionalEquations)} + * method. + *

+ * @see ExpandableFirstOrderDifferentialEquations + * @version $Id$ + * @since 3.0 + */ +public interface AdditionalEquations { + + /** Get the dimension of the additional state parameters. + * @return dimension of the additional state parameters + */ + int getDimension(); + + /** Compute the derivatives related to the additional state parameters. + * @param t current value of the independent time 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 + */ + void computeDerivatives(double t, double[] y, double[] yDot, double[] z, double[] zDot); + +} diff --git a/src/main/java/org/apache/commons/math/ode/AdditionalStateAndEquations.java b/src/main/java/org/apache/commons/math/ode/AdditionalStateAndEquations.java new file mode 100644 index 000000000..b1ba580fd --- /dev/null +++ b/src/main/java/org/apache/commons/math/ode/AdditionalStateAndEquations.java @@ -0,0 +1,88 @@ +/* + * 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. + *

+ * 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. + *

+ * + * @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. + *

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. + *

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]; + } + +} diff --git a/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderDifferentialEquations.java b/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderDifferentialEquations.java new file mode 100644 index 000000000..a89f3813e --- /dev/null +++ b/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderDifferentialEquations.java @@ -0,0 +1,272 @@ +/* + * 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. + *

+ * 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. + *

+ *

+ * 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 only the main set + * to estimate the errors and hence the step sizes. It should not 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. + *

+ *

+ * We consider that the main set always corresponds to the first equations and + * the expansion sets to the last equations. + *

+ * + * @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 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(); + } + + /** Return the dimension of the whole state vector. + *

+ * The whole state vector results in the assembly of the main set of + * equations and, if there are some, the added sets of equations. + *

+ * @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. + *

+ * 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. + *

+ * @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 time 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. + *

+ * The total current state computed by the integrator + * is dispatched here to the various additional states. + *

+ * @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); + } + +} diff --git a/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderIntegrator.java b/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderIntegrator.java new file mode 100644 index 000000000..aff61bdf8 --- /dev/null +++ b/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderIntegrator.java @@ -0,0 +1,65 @@ +/* + * 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. + *

+ * 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. + *

+ * + * @see ExpandableFirstOrderDifferentialEquations + * + * @version $Id$ + * @since 3.0 + */ + +public interface ExpandableFirstOrderIntegrator extends FirstOrderIntegrator { + + /** Integrate a set of differential equations up to the given time. + *

This method solves an Initial Value Problem (IVP).

+ *

The set of differential equations is composed of a main set, which + * can be extended by some sets of additional equations.

+ *

Since this method stores some internal state variables made + * available in its public interface during integration ({@link + * #getCurrentSignedStepsize()}), it is not thread-safe.

+ * @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 t0 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; + +} diff --git a/src/main/java/org/apache/commons/math/ode/JacobianMatrices.java b/src/main/java/org/apache/commons/math/ode/JacobianMatrices.java new file mode 100644 index 000000000..a33d2b941 --- /dev/null +++ b/src/main/java/org/apache/commons/math/ode/JacobianMatrices.java @@ -0,0 +1,471 @@ +/* + * 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.lang.reflect.Array; +import java.util.ArrayList; +import java.util.List; + +import org.apache.commons.math.exception.DimensionMismatchException; +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. + *

+ * It is intended to be packed into an {@link ExpandableFirstOrderDifferentialEquations} + * in conjunction with a main set of ODE, which may be: + *

    + *
  • a {@link FirstOrderDifferentialEquations}
    • + *
  • a {@link MainStateJacobianProvider}
    • + *
+ * In order to compute jacobian matrices with respect to some parameters of the + * main ODE set, the following parameter jacobian providers may be set: + *
    + *
  • a {@link ParameterJacobianProvider}
    • + *
  • a {@link ParameterizedODE}
    • + *
+ *

+ * + * @see ExpandableFirstOrderDifferentialEquations + * @see FirstOrderDifferentialEquations + * @see MainStateJacobianProvider + * @see ParameterJacobianProvider + * @see ParameterizedODE + * + * @version $Id$ + * @since 3.0 + */ +public class JacobianMatrices implements AdditionalEquations { + + /** Expandable first order differential equation. */ + private ExpandableFirstOrderDifferentialEquations efode; + + /** FODE without exact main jacobian computation skill. */ + private FirstOrderDifferentialEquations fode = null; + + /** FODE with exact main jacobian computation skill. */ + private MainStateJacobianProvider jode = null; + + /** FODE without exact parameter jacobian computation skill. */ + private ParameterizedODE pode = null; + + /** FODE with exact parameter jacobian computation skill. */ + private List pjp = new ArrayList();; + + /** List of parameters selected for parameter jacobian computation. */ + private List selectedParameters = null; + + /** Main state vector dimension. */ + private int stateDim; + + /** 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; + + /** Boolean for selected parameters consistency. */ + private boolean dirtyParameter = false; + + /** Simple constructor for an additional equations set computing jacobian matrices. + *

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()} + */ + public JacobianMatrices(final ExpandableFirstOrderDifferentialEquations extended, + final MainStateJacobianProvider jode) + throws IllegalArgumentException { + + checkCompatibility(extended, jode); + + efode = extended; + stateDim = efode.getMainSetDimension(); + mainJacobianInARow = new double[stateDim * stateDim]; + this.jode = jode; + efode.addAdditionalEquations(this); + setInitialMainStateJacobian(); + } + + /** Simple constructor for an additional equations set computing jacobian matrices. + *

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()} + */ + public JacobianMatrices(final ExpandableFirstOrderDifferentialEquations extended, + final FirstOrderDifferentialEquations fode) + throws IllegalArgumentException { + + checkCompatibility(extended, fode); + + 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 + */ + public void setParameterJacobianProvider(final ParameterJacobianProvider pjp) { + this.pjp.add(pjp); + } + + /** Add a parameter jacobian provider. + * @param pjp the parameterized ODE to compute by finite difference the parameter jacobian matrix + */ + public void setParameterizedODE(final ParameterizedODE pode) { + this.pode = pode; + dirtyParameter = true; + } + + /** Select the parameters to consider for jacobian matrices processing. + *

+ * Parameters must belong to the supported ones given by {@link + * Parameterizable#getParametersNames()}, so the main set of differential + * equations must be {@link Parameterizable}. + *

+ *

Note that each selection clears the previous selected parameters.

+ * + * @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(); + 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. + *

+ * Needed if and only if the main ODE set is a {@link ParameterizedODE} + * and the parameter has been {@link #selectParameters(String ...) selected} + *

+ *

+ * For pval, a non zero value of the parameter, pval * Math.sqrt(MathUtils.EPSILON) + * is a reasonable value for such a step. + *

+ *

+ * A zero value for such a step doesn't enable to compute the parameter jacobian matrix. + *

+ * @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; + } + } + if (!found) { + 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. + *

+ * Needed if and only if the main set is a {@link FirstOrderDifferentialEquations}. + *

+ *

+ * Zero values for those steps don't enable to compute the main state jacobian matrix. + *

+ * @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 + * @exception IllegalArgumentException if matrix dimensions are incorrect + */ + public void setInitialMainStateJacobian(final double[][] dYdY0) { + + // Check dimensions + checkDimension(stateDim, dYdY0); + checkDimension(stateDim, dYdY0[0]); + + // store the matrix in row major order as a single dimension array + int index = 0; + for (final double[] row : dYdY0) { + System.arraycopy(row, 0, mainJacobianInARow, index, stateDim); + index += stateDim; + } + // set initial additional state value in expandable fode + efode.setInitialAdditionalState(mainJacobianInARow, this); + } + + /** Set the initial value of the jacobian matrix with respect to one parameter. + *

The parameter must be {@link #selectParameters(String...) selected}.

+ * @param pName parameter name + * @param dYdP initial jacobian matrix w.r.t. the parameter + * @exception IllegalArgumentException if matrix dimensions are incorrect + */ + public void setInitialParameterJacobian(final String pName, final double[] dYdP) { + + // Check dimensions + checkDimension(stateDim, dYdP); + + // store the matrix in a global single dimension array + boolean found = false; + int index = 0; + 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; + } + index += stateDim; + } + if (! found) { + throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER, + pName); + } + } + + /** Set the default initial value of the jacobian matrix with respect to state. + *

dYdY0 is set to the identity matrix.

+ */ + 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. + *

The parameter must be {@link #selectParameters(String...) selected}.

+ *

dYdP is set to the null matrix.

+ * @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. + *

dYdY0 is set to the identity matrix and all dYdP are set to zero.

+ */ + 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. + */ + public void getCurrentMainSetJacobian(final double[][] dYdY0) { + + // get current state for this set of equations from the expandable fode + double[] p = efode.getCurrentAdditionalState(this); + + int index = 0; + for (int i = 0; i < stateDim; i++) { + System.arraycopy(p, index, dYdY0[i], 0, stateDim); + index += stateDim; + } + + } + + /** 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); + + int index = stateDim * stateDim; + for (ParameterConfiguration param: selectedParameters) { + if (param.getParameterName().equals(pName)) { + System.arraycopy(p, index, dYdP, 0, stateDim); + break; + } + index += stateDim; + } + + } + + /** {@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) + * @throws DimensionMismatchException if the array dimension does not match the expected one + */ + private void checkDimension(final int expected, final Object array) + throws DimensionMismatchException { + int arrayDimension = (array == null) ? 0 : Array.getLength(array); + if (arrayDimension != expected) { + throw new DimensionMismatchException(arrayDimension, expected); + } + } + +} + diff --git a/src/main/java/org/apache/commons/math/ode/MainStateJacobianProvider.java b/src/main/java/org/apache/commons/math/ode/MainStateJacobianProvider.java new file mode 100644 index 000000000..935f3685d --- /dev/null +++ b/src/main/java/org/apache/commons/math/ode/MainStateJacobianProvider.java @@ -0,0 +1,36 @@ +/* + * 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; + +/** Interface expanding {@link FirstOrderDifferentialEquations first order + * differential equations} in order to compute exactly the main state jacobian + * matrix for {@link JacobianMatrices partial derivatives equations}. + * + * @version $Id$ + * @since 3.0 + */ +public interface MainStateJacobianProvider extends FirstOrderDifferentialEquations { + + /** Compute the jacobian matrix of ODE with respect to main state. + * @param t current value of the independent time 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 dFdY placeholder array where to put the jacobian matrix of the ODE w.r.t. the main state vector + */ + void computeMainStateJacobian(double t, double[] y, double[] yDot, double[][] dFdY); + +} diff --git a/src/main/java/org/apache/commons/math/ode/MainStateJacobianWrapper.java b/src/main/java/org/apache/commons/math/ode/MainStateJacobianWrapper.java new file mode 100644 index 000000000..edb4c9670 --- /dev/null +++ b/src/main/java/org/apache/commons/math/ode/MainStateJacobianWrapper.java @@ -0,0 +1,72 @@ +/* + * 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; + } + } + +} diff --git a/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java b/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java index 3f742bdfc..0a390bc21 100644 --- a/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java +++ b/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java @@ -403,7 +403,7 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator { } /** Wrapper for differential equations, ensuring start evaluations are counted. */ - private class CountingDifferentialEquations implements ExtendedFirstOrderDifferentialEquations { + private class CountingDifferentialEquations implements FirstOrderDifferentialEquations { /** Dimension of the problem. */ private final int dimension; @@ -425,10 +425,6 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator { return dimension; } - /** {@inheritDoc} */ - public int getMainSetDimension() { - return mainSetDimension; - } } } diff --git a/src/main/java/org/apache/commons/math/ode/ParameterConfiguration.java b/src/main/java/org/apache/commons/math/ode/ParameterConfiguration.java new file mode 100644 index 000000000..a86ada92f --- /dev/null +++ b/src/main/java/org/apache/commons/math/ode/ParameterConfiguration.java @@ -0,0 +1,67 @@ +/* + * 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; + +/** Simple container pairing a parameter name with a step in order to compute + * the associated jacobian matrix by finite difference. + * + * @version $Id$ + * @since 3.0 + */ +class ParameterConfiguration implements Serializable { + + /** Serializable UID. */ + private static final long serialVersionUID = 2247518849090889379L; + + /** Parameter name. */ + private String parameterName; + + /** Parameter step for finite difference computation. */ + private double hP; + + /** Parameter name and step pair constructor. + * @param parameterName parameter name + * @param hP parameter step */ + public ParameterConfiguration(final String parameterName, final double hP) { + this.parameterName = parameterName; + this.hP = hP; + } + + /** Get parameter name. + * @return parameterName parameter name + */ + public String getParameterName() { + return parameterName; + } + + /** Get parameter step. + * @return hP parameter step + */ + public double getHP() { + return hP; + } + + /** Set parameter step. + * @param hP parameter step + */ + public void setHP(final double hP) { + this.hP = hP; + } + +} diff --git a/src/main/java/org/apache/commons/math/ode/ParameterJacobianProvider.java b/src/main/java/org/apache/commons/math/ode/ParameterJacobianProvider.java new file mode 100644 index 000000000..e9b914f37 --- /dev/null +++ b/src/main/java/org/apache/commons/math/ode/ParameterJacobianProvider.java @@ -0,0 +1,43 @@ +/* + * 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; + +/** Interface to compute exactly jacobian matrix for some parameter + * when computing {@link JacobianMatrices partial derivatives equations}. + * + * @version $Id$ + * @since 3.0 + */ +public interface ParameterJacobianProvider extends Parameterizable { + + /** Compute the jacobian matrix of ODE with respect to one parameter. + *

The parameter must be one given by {@link #getParametersNames()}.

+ * @param t current value of the independent time 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 + * 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; +} diff --git a/src/main/java/org/apache/commons/math/ode/ParameterJacobianWrapper.java b/src/main/java/org/apache/commons/math/ode/ParameterJacobianWrapper.java new file mode 100644 index 000000000..15b62b20b --- /dev/null +++ b/src/main/java/org/apache/commons/math/ode/ParameterJacobianWrapper.java @@ -0,0 +1,91 @@ +/* + * 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.Collection; +import java.util.HashMap; +import java.util.Map; + +/** Wrapper class to compute jacobian matrices by finite differences for ODE + * which do not compute them by themselves. + * + * @version $Id$ + * @since 3.0 + */ +class ParameterJacobianWrapper implements ParameterJacobianProvider { + + /** Main ODE set. */ + private final FirstOrderDifferentialEquations fode; + + /** 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. */ + private final Map 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 + * @see JacobianMatrices#setParameterStep(String, double) + */ + public ParameterJacobianWrapper(final FirstOrderDifferentialEquations fode, + final ParameterizedODE pode, + final Collection paramsAndSteps) { + this.fode = fode; + this.pode = pode; + this.hParam = new HashMap(); + + // set up parameters for jacobian computation + for (final ParameterConfiguration param : paramsAndSteps) { + final String name = param.getParameterName(); + if (pode.isSupported(name)) { + hParam.put(name, param.getHP()); + } + } + } + + /** {@inheritDoc} */ + public Collection getParametersNames() { + return pode.getParametersNames(); + } + + /** {@inheritDoc} */ + public boolean isSupported(String name) { + return pode.isSupported(name); + } + + /** {@inheritDoc} */ + public void computeParameterJacobian(double t, double[] y, double[] yDot, + String paramName, double[] dFdP) { + + final int n = fode.getDimension(); + final double[] tmpDot = new double[n]; + + // compute the jacobian df/dp w.r.t. parameter + final double p = pode.getParameter(paramName); + final double hP = hParam.get(paramName); + pode.setParameter(paramName, p + hP); + fode.computeDerivatives(t, y, tmpDot); + for (int i = 0; i < n; ++i) { + dFdP[i] = (tmpDot[i] - yDot[i]) / hP; + } + pode.setParameter(paramName, p); + + } + +} diff --git a/src/main/java/org/apache/commons/math/ode/Parameterizable.java b/src/main/java/org/apache/commons/math/ode/Parameterizable.java new file mode 100644 index 000000000..f3dad46b5 --- /dev/null +++ b/src/main/java/org/apache/commons/math/ode/Parameterizable.java @@ -0,0 +1,43 @@ +/* + * 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.Collection; + +/** This interface enables to process any parameterizable object. + * + * @version $Id$ + * @since 3.0 + */ + +public interface Parameterizable { + + /** Get the names of the supported parameters. + * @return parameters names + * @see #isSupported(String) + */ + Collection getParametersNames(); + + /** Check if a parameter is supported. + *

Supported parameters are those listed by {@link #getParametersNames()}.

+ * @param name parameter name to check + * @return true if the parameter is supported + * @see #getParametersNames() + */ + boolean isSupported(String name); + +} diff --git a/src/main/java/org/apache/commons/math/ode/ParameterizedODE.java b/src/main/java/org/apache/commons/math/ode/ParameterizedODE.java new file mode 100644 index 000000000..053d27d17 --- /dev/null +++ b/src/main/java/org/apache/commons/math/ode/ParameterizedODE.java @@ -0,0 +1,42 @@ +/* + * 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; + +/** Interface to compute by finite difference jacobian matrix for some parameter + * when computing {@link JacobianMatrices partial derivatives equations}. + * + * @version $Id$ + * @since 3.0 + */ + +public interface ParameterizedODE extends Parameterizable { + + /** Get parameter value from its name. + * @param name parameter name + * @return parameter value + * @exception IllegalArgumentException if parameter is not supported + */ + double getParameter(String name) throws IllegalArgumentException; + + /** Set the value for a given parameter. + * @param name parameter name + * @param value parameter value + * @exception IllegalArgumentException if parameter is not supported + */ + void setParameter(String name, double value) throws IllegalArgumentException; + +} diff --git a/src/main/java/org/apache/commons/math/ode/ParameterizedWrapper.java b/src/main/java/org/apache/commons/math/ode/ParameterizedWrapper.java new file mode 100644 index 000000000..6ab3efb59 --- /dev/null +++ b/src/main/java/org/apache/commons/math/ode/ParameterizedWrapper.java @@ -0,0 +1,77 @@ +/* + * 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.Collection; + +import org.apache.commons.math.exception.MathIllegalArgumentException; +import org.apache.commons.math.exception.util.LocalizedFormats; + + +/** Wrapper class enabling {@link FirstOrderDifferentialEquations basic simple} + * ODE instances to be used when processing {@link JacobianMatrices}. + * + * @version $Id$ + * @since 3.0 + */ +class ParameterizedWrapper implements ParameterizedODE { + + /** Basic FODE without parameter. */ + private final FirstOrderDifferentialEquations fode; + + /** Simple constructor. + * @param ode original first order differential equations + */ + public ParameterizedWrapper(final FirstOrderDifferentialEquations ode) { + this.fode = ode; + } + + /** {@inheritDoc} */ + public int getDimension() { + return fode.getDimension(); + } + + /** {@inheritDoc} */ + public void computeDerivatives(double t, double[] y, double[] yDot) { + fode.computeDerivatives(t, y, yDot); + } + + /** {@inheritDoc} */ + public Collection getParametersNames() { + return new ArrayList(); + } + + /** {@inheritDoc} */ + public boolean isSupported(String name) { + return false; + } + + /** {@inheritDoc} */ + public double getParameter(String name) + throws MathIllegalArgumentException { + if (!isSupported(name)) { + throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER, name); + } + return Double.NaN; + } + + /** {@inheritDoc} */ + public void setParameter(String name, double value) { + } + +} diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java index 677205902..23b4d5dfe 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java @@ -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.FirstOrderDifferentialEquations; +import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations; import org.apache.commons.math.ode.sampling.NordsieckStepInterpolator; import org.apache.commons.math.ode.sampling.StepHandler; import org.apache.commons.math.util.FastMath; @@ -189,22 +189,27 @@ public class AdamsBashforthIntegrator extends AdamsIntegrator { /** {@inheritDoc} */ @Override - public double integrate(final FirstOrderDifferentialEquations equations, - final double t0, final double[] y0, - final double t, final double[] y) + public double integrate(final ExpandableFirstOrderDifferentialEquations equations, + final double t0, final double[] z0, + final double t, final double[] z) throws MathIllegalStateException, MathIllegalArgumentException { - final int n = y0.length; - sanityChecks(equations, t0, y0, t, y); + sanityChecks(equations, t0, z0, t, z); setEquations(equations); resetEvaluations(); final boolean forward = t > t0; // 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, n); + System.arraycopy(y0, 0, y, 0, totalDim); } - final double[] yDot = new double[n]; + final double[] yDot = new double[totalDim]; // set up an interpolator sharing the integrator arrays final NordsieckStepInterpolator interpolator = new NordsieckStepInterpolator(); @@ -312,6 +317,10 @@ 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); + final double stopTime = stepStart; resetInternalState(); return stopTime; diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java index 5cae0e90d..941d6ccd4 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java @@ -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.FirstOrderDifferentialEquations; +import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations; import org.apache.commons.math.ode.MultistepIntegrator; @@ -86,7 +86,7 @@ public abstract class AdamsIntegrator extends MultistepIntegrator { /** {@inheritDoc} */ @Override - public abstract double integrate(final FirstOrderDifferentialEquations equations, + public abstract double integrate(final ExpandableFirstOrderDifferentialEquations equations, final double t0, final double[] y0, final double t, final double[] y) throws MathIllegalStateException, MathIllegalArgumentException; diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java index 7199fd6b7..42569eb18 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java @@ -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.FirstOrderDifferentialEquations; +import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations; import org.apache.commons.math.ode.sampling.NordsieckStepInterpolator; import org.apache.commons.math.ode.sampling.StepHandler; import org.apache.commons.math.util.FastMath; @@ -206,24 +206,29 @@ public class AdamsMoultonIntegrator extends AdamsIntegrator { /** {@inheritDoc} */ @Override - public double integrate(final FirstOrderDifferentialEquations equations, - final double t0, final double[] y0, - final double t, final double[] y) + public double integrate(final ExpandableFirstOrderDifferentialEquations equations, + final double t0, final double[] z0, + final double t, final double[] z) throws MathIllegalStateException, MathIllegalArgumentException { - final int n = y0.length; - sanityChecks(equations, t0, y0, t, y); + sanityChecks(equations, t0, z0, t, z); setEquations(equations); resetEvaluations(); final boolean forward = t > t0; // 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, n); + System.arraycopy(y0, 0, y, 0, totalDim); } - final double[] yDot = new double[y0.length]; - final double[] yTmp = new double[y0.length]; - final double[] predictedScaled = new double[y0.length]; + final double[] yDot = new double[totalDim]; + final double[] yTmp = new double[totalDim]; + final double[] predictedScaled = new double[totalDim]; Array2DRowRealMatrix nordsieckTmp = null; // set up two interpolators sharing the integrator arrays @@ -290,7 +295,7 @@ public class AdamsMoultonIntegrator extends AdamsIntegrator { updateHighOrderDerivativesPhase2(predictedScaled, correctedScaled, nordsieckTmp); // discrete events handling - System.arraycopy(yTmp, 0, y, 0, n); + System.arraycopy(yTmp, 0, y, 0, totalDim); interpolator.reinitialize(stepEnd, stepSize, correctedScaled, nordsieckTmp); interpolator.storeTime(stepStart); interpolator.shift(); @@ -330,6 +335,10 @@ 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); + final double stopTime = stepStart; stepStart = Double.NaN; stepSize = Double.NaN; diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java index 77ab92c0f..33f2c3e23 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java @@ -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.ExtendedFirstOrderDifferentialEquations; +import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations; import org.apache.commons.math.ode.FirstOrderDifferentialEquations; import org.apache.commons.math.util.FastMath; @@ -42,12 +42,11 @@ import org.apache.commons.math.util.FastMath; * component. The user can also use only two scalar values absTol and * relTol which will be used for all components. *

- * - *

If the Ordinary Differential Equations is an {@link ExtendedFirstOrderDifferentialEquations - * extended ODE} rather than a {@link FirstOrderDifferentialEquations basic ODE}, - * then only the {@link ExtendedFirstOrderDifferentialEquations#getMainSetDimension() - * main set} part of the state vector is used for stepsize control, not the complete - * state vector. + *

+ * If the Ordinary Differential Equations is an {@link ExpandableFirstOrderDifferentialEquations + * extended ODE} rather than a {@link FirstOrderDifferentialEquations basic ODE}, then + * only the {@link ExpandableFirstOrderDifferentialEquations#getMainSet() main part} + * of the state vector is used for stepsize control, not the complete state vector. *

* *

If the estimated error for ym+1 is such that @@ -224,18 +223,14 @@ public abstract class AdaptiveStepsizeIntegrator * @exception NumberIsTooSmallException if integration span is too small */ @Override - protected void sanityChecks(final FirstOrderDifferentialEquations equations, + protected void sanityChecks(final ExpandableFirstOrderDifferentialEquations equations, final double t0, final double[] y0, final double t, final double[] y) throws DimensionMismatchException, NumberIsTooSmallException { super.sanityChecks(equations, t0, y0, t, y); - if (equations instanceof ExtendedFirstOrderDifferentialEquations) { - mainSetDimension = ((ExtendedFirstOrderDifferentialEquations) equations).getMainSetDimension(); - } else { - mainSetDimension = equations.getDimension(); - } + mainSetDimension = equations.getMainSetDimension(); if ((vecAbsoluteTolerance != null) && (vecAbsoluteTolerance.length != mainSetDimension)) { throw new DimensionMismatchException(mainSetDimension, vecAbsoluteTolerance.length); @@ -248,7 +243,6 @@ public abstract class AdaptiveStepsizeIntegrator } /** Initialize the integration step. - * @param equations differential equations set * @param forward forward integration indicator * @param order order of the method * @param scale scaling vector for the state vector (can be shorter than state vector) @@ -259,8 +253,7 @@ public abstract class AdaptiveStepsizeIntegrator * @param yDot1 work array for the first time derivative of y1 * @return first integration step */ - public double initializeStep(final FirstOrderDifferentialEquations equations, - final boolean forward, final int order, final double[] scale, + public double initializeStep(final boolean forward, final int order, final double[] scale, final double t0, final double[] y0, final double[] yDot0, final double[] y1, final double[] yDot1) { @@ -356,7 +349,7 @@ public abstract class AdaptiveStepsizeIntegrator } /** {@inheritDoc} */ - public abstract double integrate (FirstOrderDifferentialEquations equations, + public abstract double integrate (ExpandableFirstOrderDifferentialEquations equations, double t0, double[] y0, double t, double[] y) throws MathIllegalStateException, MathIllegalArgumentException; diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java index 538f399a8..90b9d59ba 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java @@ -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.FirstOrderDifferentialEquations; +import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations; import org.apache.commons.math.ode.sampling.StepHandler; import org.apache.commons.math.util.FastMath; @@ -189,24 +189,30 @@ public abstract class EmbeddedRungeKuttaIntegrator /** {@inheritDoc} */ @Override - public double integrate(final FirstOrderDifferentialEquations equations, - final double t0, final double[] y0, - final double t, final double[] y) + public double integrate(final ExpandableFirstOrderDifferentialEquations equations, + final double t0, final double[] z0, + final double t, final double[] z) throws MathIllegalStateException, MathIllegalArgumentException { - sanityChecks(equations, t0, y0, t, y); + sanityChecks(equations, t0, z0, t, z); setEquations(equations); resetEvaluations(); final boolean forward = t > t0; // 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); + System.arraycopy(y0, 0, y, 0, totalDim); } - final double[][] yDotK = new double[stages][y0.length]; - final double[] yTmp = new double[y0.length]; - final double[] yDotTmp = new double[y0.length]; + final double[][] yDotK = new double[stages][totalDim]; + final double[] yTmp = new double[totalDim]; + final double[] yDotTmp = new double[totalDim]; // set up an interpolator sharing the integrator arrays final RungeKuttaStepInterpolator interpolator = (RungeKuttaStepInterpolator) prototype.copy(); @@ -248,7 +254,7 @@ public abstract class EmbeddedRungeKuttaIntegrator scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * FastMath.abs(y[i]); } } - hNew = initializeStep(equations, forward, getOrder(), scale, + hNew = initializeStep(forward, getOrder(), scale, stepStart, y, yDotK[0], yTmp, yDotK[1]); firstTime = false; } @@ -325,6 +331,10 @@ public abstract class EmbeddedRungeKuttaIntegrator } while (!isLastStep); + // dispatch result between main and additional states + System.arraycopy(y, 0, z, 0, z.length); + equations.setCurrentAdditionalState(y); + final double stopTime = stepStart; resetInternalState(); return stopTime; diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java index 855564ad8..ea38e62b0 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java @@ -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.FirstOrderDifferentialEquations; +import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations; 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,26 +541,32 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator { /** {@inheritDoc} */ @Override - public double integrate(final FirstOrderDifferentialEquations equations, - final double t0, final double[] y0, final double t, final double[] y) + public double integrate(final ExpandableFirstOrderDifferentialEquations equations, + final double t0, final double[] z0, final double t, final double[] z) throws MathIllegalStateException, MathIllegalArgumentException { - sanityChecks(equations, t0, y0, t, y); + sanityChecks(equations, t0, z0, t, z); setEquations(equations); resetEvaluations(); final boolean forward = t > t0; // create some internal working arrays - final double[] yDot0 = new double[y0.length]; - final double[] y1 = new double[y0.length]; - final double[] yTmp = new double[y0.length]; - final double[] yTmpDot = new double[y0.length]; + 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[][] 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[y0.length]; - y1Diag[k] = new double[y0.length]; + diagonal[k] = new double[totalDim]; + y1Diag[k] = new double[totalDim]; } final double[][][] fk = new double[sequence.length][][]; @@ -631,8 +637,7 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator { } if (firstTime) { - hNew = initializeStep(equations, forward, - 2 * targetIter + 1, scale, + hNew = initializeStep(forward, 2 * targetIter + 1, scale, stepStart, y, yDot0, yTmp, yTmpDot); } diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java index e1f5ad82c..5ead850ce 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java @@ -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.FirstOrderDifferentialEquations; +import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations; import org.apache.commons.math.ode.sampling.StepHandler; import org.apache.commons.math.util.FastMath; @@ -90,17 +90,23 @@ public abstract class RungeKuttaIntegrator extends AbstractIntegrator { } /** {@inheritDoc} */ - public double integrate(final FirstOrderDifferentialEquations equations, - final double t0, final double[] y0, - final double t, final double[] y) + public double integrate(final ExpandableFirstOrderDifferentialEquations equations, + final double t0, final double[] z0, + final double t, final double[] z) throws MathIllegalStateException, MathIllegalArgumentException { - sanityChecks(equations, t0, y0, t, y); + sanityChecks(equations, t0, z0, t, z); setEquations(equations); resetEvaluations(); final boolean forward = t > t0; // 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); @@ -179,6 +185,10 @@ 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); + final double stopTime = stepStart; stepStart = Double.NaN; stepSize = Double.NaN; diff --git a/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties b/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties index ace42b9ba..d225d8e36 100644 --- a/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties +++ b/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties @@ -279,6 +279,9 @@ UNABLE_TO_PERFORM_QR_DECOMPOSITION_ON_JACOBIAN = impossible de calculer la facto UNABLE_TO_SOLVE_SINGULAR_PROBLEM = r\u00e9solution impossible : probl\u00e8me singulier 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 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}" diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml index ba45d1144..146421723 100644 --- a/src/site/xdoc/changes.xml +++ b/src/site/xdoc/changes.xml @@ -52,6 +52,12 @@ The type attribute can be add,update,fix,remove. If the output is not quite correct, check for invisible trailing spaces! --> + + Completely revamped the computation of Jacobians in ODE. This computation is now + included in the mainstream class hierarchy, not in a separate package anymore, + and it allows adding other types of equations to a main ODE, not only variational + equations for Jacobians computation. + SimpleRegression implements UpdatingMultipleLinearRegression interface. diff --git a/src/test/java/org/apache/commons/math/ode/JacobianMatricesTest.java b/src/test/java/org/apache/commons/math/ode/JacobianMatricesTest.java new file mode 100644 index 000000000..a3ef2cb4a --- /dev/null +++ b/src/test/java/org/apache/commons/math/ode/JacobianMatricesTest.java @@ -0,0 +1,598 @@ +/* + * 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.ode.nonstiff.DormandPrince54Integrator; +import org.apache.commons.math.stat.descriptive.SummaryStatistics; +import org.apache.commons.math.util.FastMath; +import org.junit.Assert; +import org.junit.Test; + +public class JacobianMatricesTest { + + @Test + public void testLowAccuracyExternalDifferentiation() { + // this test does not really test FirstOrderIntegratorWithJacobians, + // it only shows that WITHOUT this class, attempting to recover + // the jacobians from external differentiation on simple integration + // results with low accuracy gives very poor results. In fact, + // the curves dy/dp = g(b) when b varies from 2.88 to 3.08 are + // essentially noise. + // This test is taken from Hairer, Norsett and Wanner book + // Solving Ordinary Differential Equations I (Nonstiff problems), + // the curves dy/dp = g(b) are in figure 6.5 + FirstOrderIntegrator integ = + new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 }); + double hP = 1.0e-12; + SummaryStatistics residualsP0 = new SummaryStatistics(); + SummaryStatistics residualsP1 = new SummaryStatistics(); + for (double b = 2.88; b < 3.08; b += 0.001) { + Brusselator brusselator = new Brusselator(b); + double[] y = { 1.3, b }; + integ.integrate(brusselator, 0, y, 20.0, y); + double[] yP = { 1.3, b + hP }; + integ.integrate(brusselator, 0, yP, 20.0, yP); + residualsP0.addValue((yP[0] - y[0]) / hP - brusselator.dYdP0()); + residualsP1.addValue((yP[1] - y[1]) / hP - brusselator.dYdP1()); + } + Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) > 500); + Assert.assertTrue(residualsP0.getStandardDeviation() > 30); + Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) > 700); + Assert.assertTrue(residualsP1.getStandardDeviation() > 40); + } + + @Test + public void testHighAccuracyExternalDifferentiation() { + FirstOrderIntegrator integ = + new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 1.0e-10, 1.0e-10 }); + double hP = 1.0e-12; + SummaryStatistics residualsP0 = new SummaryStatistics(); + SummaryStatistics residualsP1 = new SummaryStatistics(); + for (double b = 2.88; b < 3.08; b += 0.001) { + ParamBrusselator brusselator = new ParamBrusselator(b); + double[] y = { 1.3, b }; + integ.integrate(brusselator, 0, y, 20.0, y); + double[] yP = { 1.3, b + hP }; + brusselator.setParameter("b", b + hP); + integ.integrate(brusselator, 0, yP, 20.0, yP); + residualsP0.addValue((yP[0] - y[0]) / hP - brusselator.dYdP0()); + residualsP1.addValue((yP[1] - y[1]) / hP - brusselator.dYdP1()); + } + Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) > 0.02); + Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.03); + Assert.assertTrue(residualsP0.getStandardDeviation() > 0.003); + Assert.assertTrue(residualsP0.getStandardDeviation() < 0.004); + Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) > 0.04); + Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05); + Assert.assertTrue(residualsP1.getStandardDeviation() > 0.007); + Assert.assertTrue(residualsP1.getStandardDeviation() < 0.008); + } + + @Test + public void testInternalDifferentiation() { + ExpandableFirstOrderIntegrator 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; + SummaryStatistics residualsP0 = new SummaryStatistics(); + SummaryStatistics residualsP1 = new SummaryStatistics(); + for (double b = 2.88; b < 3.08; b += 0.001) { + ParamBrusselator brusselator = new ParamBrusselator(b); + brusselator.setParameter(ParamBrusselator.B, b); + 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.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 }); + integ.setMaxEvaluations(5000); + integ.integrate(efode, 0, z, 20.0, z); + jacob.getCurrentMainSetJacobian(dZdZ0); + jacob.getCurrentParameterJacobian(ParamBrusselator.B, dZdP); +// Assert.assertEquals(5000, integ.getMaxEvaluations()); +// Assert.assertTrue(integ.getEvaluations() > 1500); +// Assert.assertTrue(integ.getEvaluations() < 2100); +// Assert.assertEquals(4 * integ.getEvaluations(), integ.getEvaluations()); + residualsP0.addValue(dZdP[0] - brusselator.dYdP0()); + residualsP1.addValue(dZdP[1] - brusselator.dYdP1()); + } + Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.02); + Assert.assertTrue(residualsP0.getStandardDeviation() < 0.003); + Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05); + Assert.assertTrue(residualsP1.getStandardDeviation() < 0.01); + } + + @Test + public void testAnalyticalDifferentiation() { + ExpandableFirstOrderIntegrator integ = + new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 }); + SummaryStatistics residualsP0 = new SummaryStatistics(); + SummaryStatistics residualsP1 = new SummaryStatistics(); + for (double b = 2.88; b < 3.08; b += 0.001) { + Brusselator brusselator = new Brusselator(b); + 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); + jacob.setInitialParameterJacobian(Brusselator.B, new double[] { 0.0, 1.0 }); + integ.setMaxEvaluations(5000); + integ.integrate(efode, 0, z, 20.0, z); + jacob.getCurrentMainSetJacobian(dZdZ0); + jacob.getCurrentParameterJacobian(Brusselator.B, dZdP); +// Assert.assertEquals(5000, integ.getMaxEvaluations()); +// Assert.assertTrue(integ.getEvaluations() > 350); +// Assert.assertTrue(integ.getEvaluations() < 510); + residualsP0.addValue(dZdP[0] - brusselator.dYdP0()); + residualsP1.addValue(dZdP[1] - brusselator.dYdP1()); + } + Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.014); + Assert.assertTrue(residualsP0.getStandardDeviation() < 0.003); + Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05); + Assert.assertTrue(residualsP1.getStandardDeviation() < 0.01); + } + + @Test + public void testFinalResult() { + + ExpandableFirstOrderIntegrator 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); + 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)); + integ.setMaxEvaluations(5000); + + double t = 18 * FastMath.PI; + integ.integrate(efode, 0, y, t, y); + for (int i = 0; i < y.length; ++i) { + Assert.assertEquals(circle.exactY(t)[i], y[i], 1.0e-9); + } + + double[][] dydy0 = new double[2][2]; + jacob.getCurrentMainSetJacobian(dydy0); + for (int i = 0; i < dydy0.length; ++i) { + for (int j = 0; j < dydy0[i].length; ++j) { + Assert.assertEquals(circle.exactDyDy0(t)[i][j], dydy0[i][j], 1.0e-9); + } + } + double[] dydcx = new double[2]; + jacob.getCurrentParameterJacobian(Circle.CX, dydcx); + for (int i = 0; i < dydcx.length; ++i) { + Assert.assertEquals(circle.exactDyDcx(t)[i], dydcx[i], 1.0e-7); + } + double[] dydcy = new double[2]; + jacob.getCurrentParameterJacobian(Circle.CY, dydcy); + for (int i = 0; i < dydcy.length; ++i) { + Assert.assertEquals(circle.exactDyDcy(t)[i], dydcy[i], 1.0e-7); + } + double[] dydom = new double[2]; + jacob.getCurrentParameterJacobian(Circle.OMEGA, dydom); + for (int i = 0; i < dydom.length; ++i) { + Assert.assertEquals(circle.exactDyDom(t)[i], dydom[i], 1.0e-7); + } + } + + @Test + public void testParameterizable() { + + ExpandableFirstOrderIntegrator 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); +// pcircle.setParameter(ParameterizedCircle.CX, 1.0); +// 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); + 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)); + integ.setMaxEvaluations(50000); + + double t = 18 * FastMath.PI; + integ.integrate(efode, 0, y, t, y); + for (int i = 0; i < y.length; ++i) { + Assert.assertEquals(pcircle.exactY(t)[i], y[i], 1.0e-9); + } + + double[][] dydy0 = new double[2][2]; + jacob.getCurrentMainSetJacobian(dydy0); + for (int i = 0; i < dydy0.length; ++i) { + for (int j = 0; j < dydy0[i].length; ++j) { + Assert.assertEquals(pcircle.exactDyDy0(t)[i][j], dydy0[i][j], 5.0e-4); + } + } + + double[] dydp0 = new double[2]; + jacob.getCurrentParameterJacobian(Circle.CX, dydp0); + for (int i = 0; i < dydp0.length; ++i) { + Assert.assertEquals(pcircle.exactDyDcx(t)[i], dydp0[i], 5.0e-4); + } + + double[] dydp1 = new double[2]; + jacob.getCurrentParameterJacobian(Circle.OMEGA, dydp1); + for (int i = 0; i < dydp1.length; ++i) { + Assert.assertEquals(pcircle.exactDyDom(t)[i], dydp1[i], 1.0e-2); + } + } + + private static class Brusselator extends AbstractParameterizable + implements MainStateJacobianProvider, ParameterJacobianProvider { + + public static final String B = "b"; + + private double b; + + public Brusselator(double b) { + super(B); + this.b = b; + } + + public int getDimension() { + return 2; + } + + public void computeDerivatives(double t, double[] y, double[] yDot) { + double prod = y[0] * y[0] * y[1]; + yDot[0] = 1 + prod - (b + 1) * y[0]; + yDot[1] = b * y[0] - prod; + } + + public void computeMainStateJacobian(double t, double[] y, double[] yDot, + double[][] dFdY) { + double p = 2 * y[0] * y[1]; + double y02 = y[0] * y[0]; + dFdY[0][0] = p - (1 + b); + dFdY[0][1] = y02; + dFdY[1][0] = b - p; + dFdY[1][1] = -y02; + } + + public void computeParameterJacobian(double t, double[] y, double[] yDot, + String paramName, double[] dFdP) { + complainIfNotSupported(paramName); + dFdP[0] = -y[0]; + dFdP[1] = y[0]; + } + + public double dYdP0() { + return -1088.232716447743 + (1050.775747149553 + (-339.012934631828 + 36.52917025056327 * b) * b) * b; + } + + public double dYdP1() { + return 1502.824469929139 + (-1438.6974831849952 + (460.959476642384 - 49.43847385647082 * b) * b) * b; + } + + } + + private static class ParamBrusselator extends AbstractParameterizable + implements FirstOrderDifferentialEquations, ParameterizedODE { + + public static final String B = "b"; + + private double b; + + public ParamBrusselator(double b) { + super(B); + this.b = b; + } + + public int getDimension() { + return 2; + } + + /** {@inheritDoc} */ + public double getParameter(final String name) + throws IllegalArgumentException { + complainIfNotSupported(name); + return b; + } + + /** {@inheritDoc} */ + public void setParameter(final String name, final double value) + throws IllegalArgumentException { + complainIfNotSupported(name); + b = value; + } + + public void computeDerivatives(double t, double[] y, double[] yDot) { + double prod = y[0] * y[0] * y[1]; + yDot[0] = 1 + prod - (b + 1) * y[0]; + yDot[1] = b * y[0] - prod; + } + + public double dYdP0() { + return -1088.232716447743 + (1050.775747149553 + (-339.012934631828 + 36.52917025056327 * b) * b) * b; + } + + public double dYdP1() { + return 1502.824469929139 + (-1438.6974831849952 + (460.959476642384 - 49.43847385647082 * b) * b) * b; + } + + } + + /** ODE representing a point moving on a circle with provided center and angular rate. */ + private static class Circle extends AbstractParameterizable + implements MainStateJacobianProvider, ParameterJacobianProvider { + + public static final String CX = "cx"; + public static final String CY = "cy"; + public static final String OMEGA = "omega"; + + private final double[] y0; + private double cx; + private double cy; + private double omega; + + public Circle(double[] y0, double cx, double cy, double omega) { + super(CX,CY,OMEGA); + this.y0 = y0.clone(); + this.cx = cx; + this.cy = cy; + this.omega = omega; + } + + public int getDimension() { + return 2; + } + + public void computeDerivatives(double t, double[] y, double[] yDot) { + yDot[0] = omega * (cy - y[1]); + yDot[1] = omega * (y[0] - cx); + } + + public void computeMainStateJacobian(double t, double[] y, + double[] yDot, double[][] dFdY) { + dFdY[0][0] = 0; + dFdY[0][1] = -omega; + dFdY[1][0] = omega; + dFdY[1][1] = 0; + } + + public void computeParameterJacobian(double t, double[] y, double[] yDot, + String paramName, double[] dFdP) { + complainIfNotSupported(paramName); + if (paramName.equals(CX)) { + dFdP[0] = 0; + dFdP[1] = -omega; + } else if (paramName.equals(CY)) { + dFdP[0] = omega; + dFdP[1] = 0; + } else { + dFdP[0] = cy - y[1]; + dFdP[1] = y[0] - cx; + } + } + + public double[] exactY(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + double dx0 = y0[0] - cx; + double dy0 = y0[1] - cy; + return new double[] { + cx + cos * dx0 - sin * dy0, + cy + sin * dx0 + cos * dy0 + }; + } + + public double[][] exactDyDy0(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + return new double[][] { + { cos, -sin }, + { sin, cos } + }; + } + + public double[] exactDyDcx(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + return new double[] {1 - cos, -sin}; + } + + public double[] exactDyDcy(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + return new double[] {sin, 1 - cos}; + } + + public double[] exactDyDom(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + double dx0 = y0[0] - cx; + double dy0 = y0[1] - cy; + return new double[] { -t * (sin * dx0 + cos * dy0) , t * (cos * dx0 - sin * dy0) }; + } + + public double[][] exactDyDp(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + double dx0 = y0[0] - cx; + double dy0 = y0[1] - cy; + return new double[][] { + { 1 - cos, sin, -t * (sin * dx0 + cos * dy0) }, + { -sin, 1 - cos, t * (cos * dx0 - sin * dy0) } + }; + } + + public double[] exactYDot(double t) { + double oCos = omega * FastMath.cos(omega * t); + double oSin = omega * FastMath.sin(omega * t); + double dx0 = y0[0] - cx; + double dy0 = y0[1] - cy; + return new double[] { + -oSin * dx0 - oCos * dy0, + oCos * dx0 - oSin * dy0 + }; + } + + public double[][] exactDyDy0Dot(double t) { + double oCos = omega * FastMath.cos(omega * t); + double oSin = omega * FastMath.sin(omega * t); + return new double[][] { + { -oSin, -oCos }, + { oCos, -oSin } + }; + } + + public double[][] exactDyDpDot(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + double oCos = omega * cos; + double oSin = omega * sin; + double dx0 = y0[0] - cx; + double dy0 = y0[1] - cy; + return new double[][] { + { oSin, oCos, -sin * dx0 - cos * dy0 - t * ( oCos * dx0 - oSin * dy0) }, + { -oCos, oSin, cos * dx0 - sin * dy0 + t * (-oSin * dx0 - oCos * dy0) } + }; + } + + } + + /** ODE representing a point moving on a circle with provided center and angular rate. */ + private static class ParameterizedCircle extends AbstractParameterizable + implements FirstOrderDifferentialEquations, ParameterJacobianProvider, ParameterizedODE { + + public static final String CX = "cx"; + public static final String CY = "cy"; + public static final String OMEGA = "omega"; + + private final double[] y0; + private double cx; + private double cy; + private double omega; + + public ParameterizedCircle(double[] y0, double cx, double cy, double omega) { + super(CX,CY,OMEGA); + this.y0 = y0.clone(); + this.cx = cx; + this.cy = cy; + this.omega = omega; + } + + public int getDimension() { + return 2; + } + + public void computeDerivatives(double t, double[] y, double[] yDot) { + yDot[0] = omega * (cy - y[1]); + yDot[1] = omega * (y[0] - cx); + } + + public void computeParameterJacobian(double t, double[] y, double[] yDot, + String paramName, double[] dFdP) + throws IllegalArgumentException { + if (paramName.equals(CX)) { + dFdP[0] = 0; + dFdP[1] = -omega; + } else { + throw new IllegalArgumentException(); + } + } + + public double getParameter(final String name) + throws IllegalArgumentException { + if (name.equals(CY)) { + return cy; + } else if (name.equals(OMEGA)) { + return omega; + } else { + throw new IllegalArgumentException(); + } + } + + public void setParameter(final String name, final double value) + throws IllegalArgumentException { + if (name.equals(CY)) { + cy = value; + } else if (name.equals(OMEGA)) { + omega = value; + } else { + throw new IllegalArgumentException(); + } + } + + public double[] exactY(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + double dx0 = y0[0] - cx; + double dy0 = y0[1] - cy; + return new double[] { + cx + cos * dx0 - sin * dy0, + cy + sin * dx0 + cos * dy0 + }; + } + + public double[][] exactDyDy0(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + return new double[][] { + { cos, -sin }, + { sin, cos } + }; + } + + public double[] exactDyDcx(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + return new double[] {1 - cos, -sin}; + } + + public double[] exactDyDcy(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + return new double[] {sin, 1 - cos}; + } + + public double[] exactDyDom(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + double dx0 = y0[0] - cx; + double dy0 = y0[1] - cy; + return new double[] { -t * (sin * dx0 + cos * dy0) , t * (cos * dx0 - sin * dy0) }; + } + + } + +} diff --git a/src/test/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java b/src/test/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java index 6646df275..69ae26a03 100644 --- a/src/test/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java +++ b/src/test/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java @@ -256,8 +256,7 @@ public class GraggBulirschStoerIntegratorTest { } @Test - public void testUnstableDerivative() - { + public void testUnstableDerivative() { final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0); FirstOrderIntegrator integ = new GraggBulirschStoerIntegrator(0.1, 10, 1.0e-12, 0.0);