revamped Jacobians computation in ODE
git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1175409 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
39eea3de3f
commit
b7c5cd5ff7
|
|
@ -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}\""),
|
||||
|
|
|
|||
|
|
@ -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<StepHandler> 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 <I>time</I> 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);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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<String> parametersNames;
|
||||
|
||||
/** Simple constructor.
|
||||
* @param names names of the supported parameters
|
||||
*/
|
||||
protected AbstractParameterizable(final String ... names) {
|
||||
parametersNames = new ArrayList<String>();
|
||||
for (final String name : names) {
|
||||
parametersNames.add(name);
|
||||
}
|
||||
}
|
||||
|
||||
/** Simple constructor.
|
||||
* @param names names of the supported parameters
|
||||
*/
|
||||
protected AbstractParameterizable(final Collection<String> names) {
|
||||
parametersNames = new ArrayList<String>();
|
||||
parametersNames.addAll(names);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public Collection<String> 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);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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.
|
||||
* <p>
|
||||
* In some cases users may need to integrate some problem-specific equations along
|
||||
* with a main set of differential equations. One example is optimal control where
|
||||
* adjoined parameters linked to the minimized hamiltonian must be integrated.
|
||||
* </p>
|
||||
* <p>
|
||||
* This interface allows users to add such equations to a main set of {@link
|
||||
* FirstOrderDifferentialEquations first order differential equations}
|
||||
* thanks to the {@link
|
||||
* ExpandableFirstOrderDifferentialEquations#addAdditionalEquations(AdditionalEquations)}
|
||||
* method.
|
||||
* </p>
|
||||
* @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 <I>time</I> variable
|
||||
* @param y array containing the current value of the main state vector
|
||||
* @param yDot array containing the derivative of the main state vector
|
||||
* @param z array containing the current value of the additional state vector
|
||||
* @param zDot placeholder array where to put the derivative of the additional state vector
|
||||
*/
|
||||
void computeDerivatives(double t, double[] y, double[] yDot, double[] z, double[] zDot);
|
||||
|
||||
}
|
||||
|
|
@ -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.
|
||||
* <p>
|
||||
* It is a container allowing the integrator to keep constant consistency between
|
||||
* additional states and the corresponding equations. It allows to set additional
|
||||
* state values, get current additional state value and derivatives by reference
|
||||
* on the associated additional equations.
|
||||
* </p>
|
||||
*
|
||||
* @see ExpandableFirstOrderDifferentialEquations
|
||||
* @see AdditionalEquations
|
||||
*
|
||||
* @version $Id$
|
||||
* @since 3.0
|
||||
*/
|
||||
class AdditionalStateAndEquations {
|
||||
|
||||
/** Additional equations set. */
|
||||
private final AdditionalEquations addEquations;
|
||||
|
||||
/** Current additional state. */
|
||||
private double[] addState;
|
||||
|
||||
/** Current additional state derivatives. */
|
||||
private double[] addStateDot;
|
||||
|
||||
/** Create a new instance based on one set of additional equations.
|
||||
* @param addEqu additional equations.
|
||||
*/
|
||||
public AdditionalStateAndEquations(final AdditionalEquations addEqu) {
|
||||
this.addEquations = addEqu;
|
||||
}
|
||||
|
||||
/** Get a reference to the current value of the additional state.
|
||||
* <p>The array returned is a true reference to the state array, so it may be
|
||||
* used to store data into it.</>
|
||||
* @return a reference current value of the additional state.
|
||||
*/
|
||||
public double[] getAdditionalState() {
|
||||
return addState;
|
||||
}
|
||||
|
||||
/** Get a reference to the current value of the additional state derivatives.
|
||||
* <p>The array returned is a true reference to the state array, so it may be
|
||||
* used to store data into it.</>
|
||||
* @return a reference current value of the additional state derivatives.
|
||||
*/
|
||||
public double[] getAdditionalStateDot() {
|
||||
return addStateDot;
|
||||
}
|
||||
|
||||
/** Get the instance of the current additional equations.
|
||||
* @return current value of the additional equations.
|
||||
*/
|
||||
public AdditionalEquations getAdditionalEquations() {
|
||||
return addEquations;
|
||||
}
|
||||
|
||||
/** Set a value to additional state.
|
||||
* @param state additional state value.
|
||||
*/
|
||||
public void setAdditionalState(final double[] state) {
|
||||
this.addState = state.clone();
|
||||
this.addStateDot = new double[state.length];
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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.
|
||||
* <p>
|
||||
* This class extends the {@link FirstOrderDifferentialEquations}. It allows to
|
||||
* identify which part of a complete set of differential equations correspond to
|
||||
* the main set and which part correspond to the expansion sets.
|
||||
* </p>
|
||||
* <p>
|
||||
* One typical use case is the computation of the jacobian matrix for some ODE.
|
||||
* The main set of equations corresponds to the raw ODE, and we add to this set
|
||||
* another bunch of equations which represent the jacobian matrix of the main
|
||||
* set. In that case, we want the integrator to use <em>only</em> the main set
|
||||
* to estimate the errors and hence the step sizes. It should <em>not</em> use
|
||||
* the additional equations in this computation.
|
||||
* The {@link ExpandableFirstOrderIntegrator integrator} will be able to know
|
||||
* where the main set ends and so where the expansion sets begin.
|
||||
* </p>
|
||||
* <p>
|
||||
* We consider that the main set always corresponds to the first equations and
|
||||
* the expansion sets to the last equations.
|
||||
* </p>
|
||||
*
|
||||
* @see FirstOrderDifferentialEquations
|
||||
* @see JacobianMatrices
|
||||
*
|
||||
* @version $Id$
|
||||
* @since 3.0
|
||||
*/
|
||||
|
||||
public class ExpandableFirstOrderDifferentialEquations implements FirstOrderDifferentialEquations {
|
||||
|
||||
/** Main set of differential equations. */
|
||||
private final FirstOrderDifferentialEquations mainSet;
|
||||
|
||||
/** Additional sets of equations and associated states. */
|
||||
private final List<AdditionalStateAndEquations> addedSets;
|
||||
|
||||
/** Create a new instance of ExpandableFirstOrderDifferentialEquations.
|
||||
* @param fode the main set of differential equations to be integrated.
|
||||
*/
|
||||
public ExpandableFirstOrderDifferentialEquations(final FirstOrderDifferentialEquations fode) {
|
||||
this.mainSet = fode;
|
||||
this.addedSets = new ArrayList<AdditionalStateAndEquations>();
|
||||
}
|
||||
|
||||
/** Return the dimension of the whole state vector.
|
||||
* <p>
|
||||
* The whole state vector results in the assembly of the main set of
|
||||
* equations and, if there are some, the added sets of equations.
|
||||
* </p>
|
||||
* @return dimension of the whole state vector
|
||||
*/
|
||||
public int getDimension()
|
||||
throws MathIllegalArgumentException {
|
||||
int dimension = this.getMainSetDimension();
|
||||
try {
|
||||
for (AdditionalStateAndEquations stateAndEqu : addedSets) {
|
||||
dimension += stateAndEqu.getAdditionalEquations().getDimension();
|
||||
}
|
||||
return dimension;
|
||||
} catch (Exception e) {
|
||||
// TODO we should not catch Exception, and we should identify the offending additional equation
|
||||
throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_ADDITIONAL_EQUATION);
|
||||
}
|
||||
}
|
||||
|
||||
/** Return the dimension of the main set of equations.
|
||||
* <p>
|
||||
* The main set of equations represents the first part of an ODE state.
|
||||
* The error estimations and adaptive step size computation should be
|
||||
* done on this first part only, not on the final part of the state
|
||||
* which represents expansion sets of equations considered as secondary.
|
||||
* </p>
|
||||
* @return dimension of the main set of equations, must be lesser than or
|
||||
* equal to the {@link #getDimension() total dimension}
|
||||
*/
|
||||
public int getMainSetDimension() {
|
||||
return mainSet.getDimension();
|
||||
}
|
||||
|
||||
/** Return the cumulated dimension of all added sets of equations.
|
||||
* @return dimension of all added sets of equations
|
||||
* @throws IllegalArgumentException if some additional equation is unknown
|
||||
*/
|
||||
public int getAddedSetsDimension()
|
||||
throws IllegalArgumentException {
|
||||
int addDim = 0;
|
||||
try {
|
||||
for (AdditionalStateAndEquations stateAndEqu : addedSets) {
|
||||
addDim += stateAndEqu.getAdditionalEquations().getDimension();
|
||||
}
|
||||
return addDim;
|
||||
} catch (Exception e) {
|
||||
// TODO we should not catch Exception, and we should identify the offending additional equation
|
||||
throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_ADDITIONAL_EQUATION);
|
||||
}
|
||||
}
|
||||
|
||||
/** Return the dimension of one added set of equations.
|
||||
* @param addEqu Additional equations used as a reference for selection
|
||||
* @return dimension of the added set of equations
|
||||
* @throws IllegalArgumentException if additional equation is unknown
|
||||
*/
|
||||
public int getAddedSetDimension(final AdditionalEquations addEqu) {
|
||||
return selectStateAndEquations(addEqu).getAdditionalEquations().getDimension();
|
||||
}
|
||||
|
||||
/** Get the current time derivative of the total state vector.
|
||||
* @param t current value of the independent <I>time</I> variable
|
||||
* @param y array containing the current value of the state vector
|
||||
* @param yDot placeholder array where to put the time derivative of the state vector
|
||||
*/
|
||||
public void computeDerivatives(final double t, final double[] y, final double[] yDot) {
|
||||
|
||||
// Add contribution for the main set of equations
|
||||
int index = getMainSetDimension();
|
||||
double[] m = new double[index];
|
||||
double[] mDot = new double[index];
|
||||
// update current main state
|
||||
System.arraycopy(y, 0, m, 0, index);
|
||||
// compute derivatives
|
||||
mainSet.computeDerivatives(t, m, mDot);
|
||||
// update main state contribution in global array
|
||||
System.arraycopy(mDot, 0, yDot, 0, index);
|
||||
|
||||
// Add contribution for additional equations
|
||||
for (final AdditionalStateAndEquations stateAndEqu : addedSets) {
|
||||
final double[] p = stateAndEqu.getAdditionalState();
|
||||
final double[] pDot = stateAndEqu.getAdditionalStateDot();
|
||||
|
||||
// update current additional state
|
||||
System.arraycopy(y, index, p, 0, p.length);
|
||||
|
||||
// compute additional derivatives
|
||||
stateAndEqu.getAdditionalEquations().computeDerivatives(t, m, mDot, p, pDot);
|
||||
|
||||
// update each additional state contribution in global array
|
||||
System.arraycopy(pDot, 0, yDot, index, p.length);
|
||||
|
||||
// incrementing index
|
||||
index += p.length;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/** Add a set of user-specified equations to be integrated along with
|
||||
* the main set of equations.
|
||||
*
|
||||
* @param addEqu additional equations
|
||||
* @see #setInitialAdditionalState(double[], AdditionalEquations)
|
||||
* @see #getCurrentAdditionalState(AdditionalEquations)
|
||||
*/
|
||||
public void addAdditionalEquations(final AdditionalEquations addEqu) {
|
||||
addedSets.add(new AdditionalStateAndEquations(addEqu));
|
||||
}
|
||||
|
||||
/** Get the instance of the main set of equations.
|
||||
* @return current value of the main set of equations.
|
||||
*/
|
||||
public FirstOrderDifferentialEquations getMainSet() {
|
||||
return mainSet;
|
||||
}
|
||||
|
||||
/** Set initial additional state.
|
||||
* @param addState additional state
|
||||
* @param addEqu additional equations used as a reference for selection
|
||||
* @throws IllegalArgumentException if additional equation is unknown
|
||||
*/
|
||||
public void setInitialAdditionalState(final double[] addState, final AdditionalEquations addEqu) {
|
||||
selectStateAndEquations(addEqu).setAdditionalState(addState);
|
||||
}
|
||||
|
||||
/** Set current additional state.
|
||||
* <p>
|
||||
* The total current state computed by the integrator
|
||||
* is dispatched here to the various additional states.
|
||||
* </p>
|
||||
* @param currentState total current state
|
||||
* @throws IllegalArgumentException if additional equation is unknown
|
||||
*/
|
||||
public void setCurrentAdditionalState(final double[] currentState)
|
||||
throws IllegalArgumentException {
|
||||
int index = getMainSetDimension();
|
||||
try {
|
||||
for (AdditionalStateAndEquations stateAndEqu : addedSets) {
|
||||
final int addDim = stateAndEqu.getAdditionalEquations().getDimension();
|
||||
final double[] addState = new double[addDim];
|
||||
System.arraycopy(currentState, index, addState, 0, addDim);
|
||||
stateAndEqu.setAdditionalState(addState);
|
||||
index += addDim;
|
||||
}
|
||||
} catch (Exception e) {
|
||||
// TODO we should not catch Exception, and we should identify the offending additional equation
|
||||
throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_ADDITIONAL_EQUATION);
|
||||
}
|
||||
}
|
||||
|
||||
/** Get current additional state.
|
||||
* @param addEqu additional equations used as a reference for selection
|
||||
* @return current additional state
|
||||
* @throws IllegalArgumentException if additional equation is unknown
|
||||
*/
|
||||
public double[] getCurrentAdditionalState(final AdditionalEquations addEqu) {
|
||||
return selectStateAndEquations(addEqu).getAdditionalState();
|
||||
}
|
||||
|
||||
/** Get all current additional states accumulated.
|
||||
* @return current additional states
|
||||
* @throws IllegalArgumentException if additional equation is unknown
|
||||
*/
|
||||
public double[] getCurrentAdditionalStates()
|
||||
throws IllegalArgumentException {
|
||||
int index = 0;
|
||||
final double[] cumulState = new double[getAddedSetsDimension()];
|
||||
try {
|
||||
for (AdditionalStateAndEquations stateAndEqu : addedSets) {
|
||||
final int addDim = stateAndEqu.getAdditionalEquations().getDimension();
|
||||
final double[] addState = stateAndEqu.getAdditionalState();
|
||||
System.arraycopy(addState, 0, cumulState, index, addDim);
|
||||
index += addDim;
|
||||
}
|
||||
return cumulState;
|
||||
} catch (Exception e) {
|
||||
// TODO we should not catch Exception, and we should identify the offending additional equation
|
||||
throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_ADDITIONAL_EQUATION);
|
||||
}
|
||||
}
|
||||
|
||||
/** Select additional state and equations pair in the list.
|
||||
* @param addEqu Additional equations used as a reference for selection
|
||||
* @return additional state and equations pair
|
||||
* @throws IllegalArgumentException if additional equation is unknown
|
||||
*/
|
||||
private AdditionalStateAndEquations selectStateAndEquations(final AdditionalEquations addEqu)
|
||||
throws IllegalArgumentException {
|
||||
for (AdditionalStateAndEquations stateAndEqu : addedSets) {
|
||||
if (stateAndEqu.getAdditionalEquations() == addEqu) {
|
||||
return stateAndEqu;
|
||||
}
|
||||
}
|
||||
// TODO we should not catch Exception, and we should identify the offending additional equation
|
||||
throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_ADDITIONAL_EQUATION);
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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.
|
||||
* <p>
|
||||
* The classes devoted to solve expandable first order differential equations
|
||||
* should implement this interface. The problems which can be handled should
|
||||
* implement the {@link ExpandableFirstOrderDifferentialEquations} interface.
|
||||
* </p>
|
||||
*
|
||||
* @see ExpandableFirstOrderDifferentialEquations
|
||||
*
|
||||
* @version $Id$
|
||||
* @since 3.0
|
||||
*/
|
||||
|
||||
public interface ExpandableFirstOrderIntegrator extends FirstOrderIntegrator {
|
||||
|
||||
/** Integrate a set of differential equations up to the given time.
|
||||
* <p>This method solves an Initial Value Problem (IVP).</p>
|
||||
* <p>The set of differential equations is composed of a main set, which
|
||||
* can be extended by some sets of additional equations.</p>
|
||||
* <p>Since this method stores some internal state variables made
|
||||
* available in its public interface during integration ({@link
|
||||
* #getCurrentSignedStepsize()}), it is <em>not</em> thread-safe.</p>
|
||||
* @param equations complete set of differential equations to integrate
|
||||
* @param t0 initial time
|
||||
* @param y0 initial value of the main state vector at t0
|
||||
* @param t target time for the integration
|
||||
* (can be set to a value smaller than <code>t0</code> for backward integration)
|
||||
* @param y placeholder where to put the main state vector at each successful
|
||||
* step (and hence at the end of integration), can be the same object as y0
|
||||
* @return stop time, will be the same as target time if integration reached its
|
||||
* target, but may be different if some {@link
|
||||
* org.apache.commons.math.ode.events.EventHandler} stops it at some point.
|
||||
* @throws MathIllegalStateException if the integrator cannot perform integration
|
||||
* @throws MathIllegalArgumentException if integration parameters are wrong (typically
|
||||
* too small integration span)
|
||||
*/
|
||||
double integrate(ExpandableFirstOrderDifferentialEquations equations,
|
||||
double t0, double[] y0, double t, double[] y)
|
||||
throws MathIllegalStateException, MathIllegalArgumentException;
|
||||
|
||||
}
|
||||
|
|
@ -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.
|
||||
* <p>
|
||||
* It is intended to be packed into an {@link ExpandableFirstOrderDifferentialEquations}
|
||||
* in conjunction with a main set of ODE, which may be:
|
||||
* <ul>
|
||||
* <li>a {@link FirstOrderDifferentialEquations}</li>
|
||||
* <li>a {@link MainStateJacobianProvider}</li>
|
||||
* </ul>
|
||||
* In order to compute jacobian matrices with respect to some parameters of the
|
||||
* main ODE set, the following parameter jacobian providers may be set:
|
||||
* <ul>
|
||||
* <li>a {@link ParameterJacobianProvider}</li>
|
||||
* <li>a {@link ParameterizedODE}</li>
|
||||
* </ul>
|
||||
* </p>
|
||||
*
|
||||
* @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<ParameterJacobianProvider> pjp = new ArrayList<ParameterJacobianProvider>();;
|
||||
|
||||
/** List of parameters selected for parameter jacobian computation. */
|
||||
private List<ParameterConfiguration> 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.
|
||||
* <p>This additional equations set is added internally to the expandable
|
||||
* first order differential equations set thanks to the
|
||||
* {@link ExpandableFirstOrderDifferentialEquations#addAdditionalEquations(AdditionalEquations)}
|
||||
* method.
|
||||
* @param extended the expandable first order differential equations set
|
||||
* @param jode the main first order differential equations set to extend
|
||||
* @exception IllegalArgumentException if jode does not match the main set to be extended given by
|
||||
* {@link ExpandableFirstOrderDifferentialEquations#getMainSet() extended.getMainSet()}
|
||||
*/
|
||||
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.
|
||||
* <p>This additional equations set is added internally to the expandable
|
||||
* first order differential equations set thanks to the
|
||||
* {@link ExpandableFirstOrderDifferentialEquations#addAdditionalEquations(AdditionalEquations)}
|
||||
* method.
|
||||
* @param extended the expandable first order differential equations set
|
||||
* @param fode the main first order differential equations set to extend
|
||||
* @exception IllegalArgumentException if fode does not match the main set to be extended given by
|
||||
* {@link ExpandableFirstOrderDifferentialEquations#getMainSet() extended.getMainSet()}
|
||||
*/
|
||||
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.
|
||||
* <p>
|
||||
* Parameters must belong to the supported ones given by {@link
|
||||
* Parameterizable#getParametersNames()}, so the main set of differential
|
||||
* equations must be {@link Parameterizable}.
|
||||
* </p>
|
||||
* <p>Note that each selection clears the previous selected parameters.</p>
|
||||
*
|
||||
* @param parameters parameters to consider for jacobian matrices processing
|
||||
* @exception IllegalArgumentException if one parameter is not supported
|
||||
*/
|
||||
public void selectParameters(final String... parameters) throws IllegalArgumentException {
|
||||
|
||||
selectedParameters = new ArrayList<ParameterConfiguration>();
|
||||
for (String param : parameters) {
|
||||
selectedParameters.add(new ParameterConfiguration(param, Double.NaN));
|
||||
}
|
||||
paramDim = parameters.length;
|
||||
parameterJacobiansInARow = new double[paramDim * stateDim];
|
||||
setInitialParameterJacobians();
|
||||
}
|
||||
|
||||
/** Set the step associated to a parameter in order to compute by finite
|
||||
* difference the jacobian matrix.
|
||||
* <p>
|
||||
* Needed if and only if the main ODE set is a {@link ParameterizedODE}
|
||||
* and the parameter has been {@link #selectParameters(String ...) selected}
|
||||
* </p>
|
||||
* <p>
|
||||
* For pval, a non zero value of the parameter, pval * Math.sqrt(MathUtils.EPSILON)
|
||||
* is a reasonable value for such a step.
|
||||
* </p>
|
||||
* <p>
|
||||
* A zero value for such a step doesn't enable to compute the parameter jacobian matrix.
|
||||
* </p>
|
||||
* @param parameter parameter to consider for jacobian processing
|
||||
* @param hP step for jacobian finite difference computation w.r.t. the specified parameter
|
||||
* @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.
|
||||
* <p>
|
||||
* Needed if and only if the main set is a {@link FirstOrderDifferentialEquations}.
|
||||
* </p>
|
||||
* <p>
|
||||
* Zero values for those steps don't enable to compute the main state jacobian matrix.
|
||||
* </p>
|
||||
* @param hY step used for finite difference computation with respect to state vector
|
||||
* @exception IllegalArgumentException if the hY has not the dimension of the main state
|
||||
* given by {@link ExpandableFirstOrderDifferentialEquations#getMainSetDimension()}
|
||||
*/
|
||||
public void setMainStateSteps(final double[] hY) {
|
||||
|
||||
if (fode != null) {
|
||||
// Check dimension
|
||||
checkDimension(stateDim, hY);
|
||||
this.hY = hY.clone();
|
||||
dirtyMainState = true;
|
||||
}
|
||||
}
|
||||
|
||||
/** Set the initial value of the jacobian matrix with respect to state.
|
||||
* @param dYdY0 initial jacobian matrix w.r.t. state
|
||||
* @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.
|
||||
* <p>The parameter must be {@link #selectParameters(String...) selected}.</p>
|
||||
* @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.
|
||||
* <p>dYdY0 is set to the identity matrix.</p>
|
||||
*/
|
||||
public void setInitialMainStateJacobian() {
|
||||
final double[][] dYdY0 = new double[stateDim][stateDim];
|
||||
for (int i = 0; i < stateDim; ++i) {
|
||||
dYdY0[i][i] = 1.0;
|
||||
}
|
||||
setInitialMainStateJacobian(dYdY0);
|
||||
}
|
||||
|
||||
/** Set the default initial value of the jacobian matrix with respect to one parameter.
|
||||
* <p>The parameter must be {@link #selectParameters(String...) selected}.</p>
|
||||
* <p>dYdP is set to the null matrix.</p>
|
||||
* @param pName parameter name
|
||||
*/
|
||||
public void setInitialParameterJacobian(final String pName) {
|
||||
setInitialParameterJacobian(pName, new double[stateDim]);
|
||||
}
|
||||
|
||||
/** Set the default initial values of jacobian matrices with respect to all parameters.
|
||||
*/
|
||||
public void setInitialParameterJacobians() {
|
||||
for (ParameterConfiguration param: selectedParameters) {
|
||||
setInitialParameterJacobian(param.getParameterName());
|
||||
}
|
||||
}
|
||||
|
||||
/** Set default initial values for jacobian matrices.
|
||||
* <p>dYdY0 is set to the identity matrix and all dYdP are set to zero.</p>
|
||||
*/
|
||||
public void setInitialJacobians() {
|
||||
setInitialMainStateJacobian();
|
||||
setInitialParameterJacobians();
|
||||
}
|
||||
|
||||
/** Get the current value of the jacobian matrix with respect to state.
|
||||
* @param dYdY0 current jacobian matrix with respect to state.
|
||||
*/
|
||||
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);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
|
@ -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 <I>time</I> variable
|
||||
* @param y array containing the current value of the main state vector
|
||||
* @param yDot array containing the current value of the time derivative of the main state vector
|
||||
* @param 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);
|
||||
|
||||
}
|
||||
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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.
|
||||
* <p>The parameter must be one given by {@link #getParametersNames()}.</p>
|
||||
* @param t current value of the independent <I>time</I> variable
|
||||
* @param y array containing the current value of the main state vector
|
||||
* @param yDot array containing the current value of the time derivative
|
||||
* of the main state vector
|
||||
* @param paramName name of the parameter to consider
|
||||
* @param dFdP placeholder array where to put the jacobian matrix of the
|
||||
* 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;
|
||||
}
|
||||
|
|
@ -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<String, Double> hParam;
|
||||
|
||||
/** Wrap a {@link ParameterizedODE} into a {@link ParameterJacobianProvider}.
|
||||
* @param fode main first order differential equations set
|
||||
* @param pode additional problem, without parametre jacobian computation skill
|
||||
* @param paramsAndSteps parameters and steps to compute the jacobians df/dp
|
||||
* @see JacobianMatrices#setParameterStep(String, double)
|
||||
*/
|
||||
public ParameterJacobianWrapper(final FirstOrderDifferentialEquations fode,
|
||||
final ParameterizedODE pode,
|
||||
final Collection<ParameterConfiguration> paramsAndSteps) {
|
||||
this.fode = fode;
|
||||
this.pode = pode;
|
||||
this.hParam = new HashMap<String, Double>();
|
||||
|
||||
// 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<String> 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);
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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<String> getParametersNames();
|
||||
|
||||
/** Check if a parameter is supported.
|
||||
* <p>Supported parameters are those listed by {@link #getParametersNames()}.</p>
|
||||
* @param name parameter name to check
|
||||
* @return true if the parameter is supported
|
||||
* @see #getParametersNames()
|
||||
*/
|
||||
boolean isSupported(String name);
|
||||
|
||||
}
|
||||
|
|
@ -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;
|
||||
|
||||
}
|
||||
|
|
@ -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<String> getParametersNames() {
|
||||
return new ArrayList<String>();
|
||||
}
|
||||
|
||||
/** {@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) {
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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.
|
||||
* </p>
|
||||
*
|
||||
* <p>If the Ordinary Differential Equations is an {@link ExtendedFirstOrderDifferentialEquations
|
||||
* extended ODE} rather than a {@link FirstOrderDifferentialEquations basic ODE},
|
||||
* then <em>only</em> the {@link ExtendedFirstOrderDifferentialEquations#getMainSetDimension()
|
||||
* main set} part of the state vector is used for stepsize control, not the complete
|
||||
* state vector.
|
||||
* <p>
|
||||
* If the Ordinary Differential Equations is an {@link ExpandableFirstOrderDifferentialEquations
|
||||
* extended ODE} rather than a {@link FirstOrderDifferentialEquations basic ODE}, then
|
||||
* <em>only</em> the {@link ExpandableFirstOrderDifferentialEquations#getMainSet() main part}
|
||||
* of the state vector is used for stepsize control, not the complete state vector.
|
||||
* </p>
|
||||
*
|
||||
* <p>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;
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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}"
|
||||
|
|
|
|||
|
|
@ -52,6 +52,12 @@ The <action> type attribute can be add,update,fix,remove.
|
|||
If the output is not quite correct, check for invisible trailing spaces!
|
||||
-->
|
||||
<release version="3.0" date="TBD" description="TBD">
|
||||
<action dev="luc" type="fix" issue="MATH-380" due-to="Pascal Parraud">
|
||||
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.
|
||||
</action>
|
||||
<action dev="gregs" type="update" issue="MATH-607">
|
||||
SimpleRegression implements UpdatingMultipleLinearRegression interface.
|
||||
</action>
|
||||
|
|
|
|||
|
|
@ -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) };
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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);
|
||||
|
|
|
|||
Loading…
Reference in New Issue