Return-Path: X-Original-To: apmail-commons-commits-archive@minotaur.apache.org Delivered-To: apmail-commons-commits-archive@minotaur.apache.org Received: from mail.apache.org (hermes.apache.org [140.211.11.3]) by minotaur.apache.org (Postfix) with SMTP id 5D2EE92BF for ; Sun, 25 Sep 2011 15:05:07 +0000 (UTC) Received: (qmail 36849 invoked by uid 500); 25 Sep 2011 15:05:07 -0000 Delivered-To: apmail-commons-commits-archive@commons.apache.org Received: (qmail 36789 invoked by uid 500); 25 Sep 2011 15:05:07 -0000 Mailing-List: contact commits-help@commons.apache.org; run by ezmlm Precedence: bulk List-Help: List-Unsubscribe: List-Post: List-Id: Reply-To: dev@commons.apache.org Delivered-To: mailing list commits@commons.apache.org Received: (qmail 36782 invoked by uid 99); 25 Sep 2011 15:05:06 -0000 Received: from athena.apache.org (HELO athena.apache.org) (140.211.11.136) by apache.org (qpsmtpd/0.29) with ESMTP; Sun, 25 Sep 2011 15:05:06 +0000 X-ASF-Spam-Status: No, hits=-2000.0 required=5.0 tests=ALL_TRUSTED X-Spam-Check-By: apache.org Received: from [140.211.11.4] (HELO eris.apache.org) (140.211.11.4) by apache.org (qpsmtpd/0.29) with ESMTP; Sun, 25 Sep 2011 15:05:02 +0000 Received: from eris.apache.org (localhost [127.0.0.1]) by eris.apache.org (Postfix) with ESMTP id 6D1A723888FD for ; Sun, 25 Sep 2011 15:04:42 +0000 (UTC) Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Subject: svn commit: r1175409 [1/2] - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/exception/util/ main/java/org/apache/commons/math/ode/ main/java/org/apache/commons/math/ode/nonstiff/ main/resources/META-INF/localization/ site/xdoc/ te... Date: Sun, 25 Sep 2011 15:04:40 -0000 To: commits@commons.apache.org From: luc@apache.org X-Mailer: svnmailer-1.0.8-patched Message-Id: <20110925150442.6D1A723888FD@eris.apache.org> Author: luc Date: Sun Sep 25 15:04:39 2011 New Revision: 1175409 URL: http://svn.apache.org/viewvc?rev=1175409&view=rev Log: revamped Jacobians computation in ODE Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractParameterizable.java (with props) commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalEquations.java (with props) commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalStateAndEquations.java (with props) commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderDifferentialEquations.java (with props) commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderIntegrator.java (with props) commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/JacobianMatrices.java (with props) commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MainStateJacobianProvider.java (with props) commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MainStateJacobianWrapper.java (with props) commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterConfiguration.java (with props) commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianProvider.java (with props) commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianWrapper.java (with props) commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/Parameterizable.java (with props) commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedODE.java (with props) commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedWrapper.java (with props) commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/JacobianMatricesTest.java (with props) Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties commons/proper/math/trunk/src/site/xdoc/changes.xml commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java?rev=1175409&r1=1175408&r2=1175409&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java Sun Sep 25 15:04:39 2011 @@ -315,6 +315,9 @@ public enum LocalizedFormats implements 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}\""), Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java?rev=1175409&r1=1175408&r2=1175409&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java Sun Sep 25 15:04:39 2011 @@ -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.Math * @version $Id$ * @since 2.0 */ -public abstract class AbstractIntegrator implements FirstOrderIntegrator { +public abstract class AbstractIntegrator implements ExpandableFirstOrderIntegrator { /** Step handler. */ protected Collection stepHandlers; @@ -76,7 +77,7 @@ public abstract class AbstractIntegrator 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 * @param equations differential equations to integrate * @see #computeDerivatives(double, double[], double[]) */ - protected void setEquations(final FirstOrderDifferentialEquations equations) { + protected void setEquations(final ExpandableFirstOrderDifferentialEquations equations) { this.equations = equations; } + /** {@inheritDoc} */ + public double integrate(FirstOrderDifferentialEquations equations, + double t0, double[] y0, double t, double[] y) + throws MathIllegalStateException, MathIllegalArgumentException { + return integrate(new ExpandableFirstOrderDifferentialEquations(equations), t0, y0, t, y); + } + /** Compute the derivatives and check the number of evaluations. * @param t current value of the independent time variable * @param y array containing the current value of the state vector @@ -336,16 +344,16 @@ public abstract class AbstractIntegrator * @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); } Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractParameterizable.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractParameterizable.java?rev=1175409&view=auto ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractParameterizable.java (added) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractParameterizable.java Sun Sep 25 15:04:39 2011 @@ -0,0 +1,81 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math.ode; + +import java.util.ArrayList; +import java.util.Collection; + +import org.apache.commons.math.exception.MathIllegalArgumentException; +import org.apache.commons.math.exception.util.LocalizedFormats; + +/** This abstract class provides boilerplate parameters list. + * + * @version $Id$ + * @since 3.0 + */ + +public abstract class AbstractParameterizable implements Parameterizable { + + /** List of the parameters names. */ + private final Collection parametersNames; + + /** Simple constructor. + * @param names names of the supported parameters + */ + protected AbstractParameterizable(final String ... names) { + parametersNames = new ArrayList(); + for (final String name : names) { + parametersNames.add(name); + } + } + + /** Simple constructor. + * @param names names of the supported parameters + */ + protected AbstractParameterizable(final Collection names) { + parametersNames = new ArrayList(); + parametersNames.addAll(names); + } + + /** {@inheritDoc} */ + public Collection getParametersNames() { + return parametersNames; + } + + /** {@inheritDoc} */ + public boolean isSupported(final String name) { + for (final String supportedName : parametersNames) { + if (supportedName.equals(name)) { + return true; + } + } + return false; + } + + /** Check if a parameter is supported and throw an IllegalArgumentException if not. + * @param name name of the parameter to check + * @exception MathIllegalArgumentException if the parameter is not supported + * @see #isSupported(String) + */ + public void complainIfNotSupported(final String name) + throws MathIllegalArgumentException { + if (!isSupported(name)) { + throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER, name); + } + } + +} Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractParameterizable.java ------------------------------------------------------------------------------ svn:eol-style = native Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractParameterizable.java ------------------------------------------------------------------------------ svn:keywords = "Author Date Id Revision" Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalEquations.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalEquations.java?rev=1175409&view=auto ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalEquations.java (added) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalEquations.java Sun Sep 25 15:04:39 2011 @@ -0,0 +1,55 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math.ode; + +/** + * This interface allows users to add their own differential equations to a main + * set of differential equations. + *

+ * In some cases users may need to integrate some problem-specific equations along + * with a main set of differential equations. One example is optimal control where + * adjoined parameters linked to the minimized hamiltonian must be integrated. + *

+ *

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

+ * @see ExpandableFirstOrderDifferentialEquations + * @version $Id$ + * @since 3.0 + */ +public interface AdditionalEquations { + + /** Get the dimension of the additional state parameters. + * @return dimension of the additional state parameters + */ + int getDimension(); + + /** Compute the derivatives related to the additional state parameters. + * @param t current value of the independent time variable + * @param y array containing the current value of the main state vector + * @param yDot array containing the derivative of the main state vector + * @param z array containing the current value of the additional state vector + * @param zDot placeholder array where to put the derivative of the additional state vector + */ + void computeDerivatives(double t, double[] y, double[] yDot, double[] z, double[] zDot); + +} Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalEquations.java ------------------------------------------------------------------------------ svn:eol-style = native Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalEquations.java ------------------------------------------------------------------------------ svn:keywords = "Author Date Id Revision" Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalStateAndEquations.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalStateAndEquations.java?rev=1175409&view=auto ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalStateAndEquations.java (added) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalStateAndEquations.java Sun Sep 25 15:04:39 2011 @@ -0,0 +1,88 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math.ode; + + +/** + * This class is a container for additional state parameters and their associated + * evolution equation. + *

+ * It is a container allowing the integrator to keep constant consistency between + * additional states and the corresponding equations. It allows to set additional + * state values, get current additional state value and derivatives by reference + * on the associated additional equations. + *

+ * + * @see ExpandableFirstOrderDifferentialEquations + * @see AdditionalEquations + * + * @version $Id$ + * @since 3.0 + */ +class AdditionalStateAndEquations { + + /** Additional equations set. */ + private final AdditionalEquations addEquations; + + /** Current additional state. */ + private double[] addState; + + /** Current additional state derivatives. */ + private double[] addStateDot; + + /** Create a new instance based on one set of additional equations. + * @param addEqu additional equations. + */ + public AdditionalStateAndEquations(final AdditionalEquations addEqu) { + this.addEquations = addEqu; + } + + /** Get a reference to the current value of the additional state. + *

The array returned is a true reference to the state array, so it may be + * used to store data into it. + * @return a reference current value of the additional state. + */ + public double[] getAdditionalState() { + return addState; + } + + /** Get a reference to the current value of the additional state derivatives. + *

The array returned is a true reference to the state array, so it may be + * used to store data into it. + * @return a reference current value of the additional state derivatives. + */ + public double[] getAdditionalStateDot() { + return addStateDot; + } + + /** Get the instance of the current additional equations. + * @return current value of the additional equations. + */ + public AdditionalEquations getAdditionalEquations() { + return addEquations; + } + + /** Set a value to additional state. + * @param state additional state value. + */ + public void setAdditionalState(final double[] state) { + this.addState = state.clone(); + this.addStateDot = new double[state.length]; + } + +} Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalStateAndEquations.java ------------------------------------------------------------------------------ svn:eol-style = native Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalStateAndEquations.java ------------------------------------------------------------------------------ svn:keywords = "Author Date Id Revision" Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderDifferentialEquations.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderDifferentialEquations.java?rev=1175409&view=auto ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderDifferentialEquations.java (added) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderDifferentialEquations.java Sun Sep 25 15:04:39 2011 @@ -0,0 +1,272 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math.ode; + +import java.util.ArrayList; +import java.util.List; + +import org.apache.commons.math.exception.MathIllegalArgumentException; +import org.apache.commons.math.exception.util.LocalizedFormats; + + +/** + * This class represents a combined set of first order differential equations, + * with at least a main set of equations expandable by some sets of additional + * equations. + *

+ * This class extends the {@link FirstOrderDifferentialEquations}. It allows to + * identify which part of a complete set of differential equations correspond to + * the main set and which part correspond to the expansion sets. + *

+ *

+ * One typical use case is the computation of the jacobian matrix for some ODE. + * The main set of equations corresponds to the raw ODE, and we add to this set + * another bunch of equations which represent the jacobian matrix of the main + * set. In that case, we want the integrator to use only the main set + * to estimate the errors and hence the step sizes. It should not use + * the additional equations in this computation. + * The {@link ExpandableFirstOrderIntegrator integrator} will be able to know + * where the main set ends and so where the expansion sets begin. + *

+ *

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

+ * + * @see FirstOrderDifferentialEquations + * @see JacobianMatrices + * + * @version $Id$ + * @since 3.0 + */ + +public class ExpandableFirstOrderDifferentialEquations implements FirstOrderDifferentialEquations { + + /** Main set of differential equations. */ + private final FirstOrderDifferentialEquations mainSet; + + /** Additional sets of equations and associated states. */ + private final List addedSets; + + /** Create a new instance of ExpandableFirstOrderDifferentialEquations. + * @param fode the main set of differential equations to be integrated. + */ + public ExpandableFirstOrderDifferentialEquations(final FirstOrderDifferentialEquations fode) { + this.mainSet = fode; + this.addedSets = new ArrayList(); + } + + /** Return the dimension of the whole state vector. + *

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

+ * @return dimension of the whole state vector + */ + public int getDimension() + throws MathIllegalArgumentException { + int dimension = this.getMainSetDimension(); + try { + for (AdditionalStateAndEquations stateAndEqu : addedSets) { + dimension += stateAndEqu.getAdditionalEquations().getDimension(); + } + return dimension; + } catch (Exception e) { + // TODO we should not catch Exception, and we should identify the offending additional equation + throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_ADDITIONAL_EQUATION); + } + } + + /** Return the dimension of the main set of equations. + *

+ * The main set of equations represents the first part of an ODE state. + * The error estimations and adaptive step size computation should be + * done on this first part only, not on the final part of the state + * which represents expansion sets of equations considered as secondary. + *

+ * @return dimension of the main set of equations, must be lesser than or + * equal to the {@link #getDimension() total dimension} + */ + public int getMainSetDimension() { + return mainSet.getDimension(); + } + + /** Return the cumulated dimension of all added sets of equations. + * @return dimension of all added sets of equations + * @throws IllegalArgumentException if some additional equation is unknown + */ + public int getAddedSetsDimension() + throws IllegalArgumentException { + int addDim = 0; + try { + for (AdditionalStateAndEquations stateAndEqu : addedSets) { + addDim += stateAndEqu.getAdditionalEquations().getDimension(); + } + return addDim; + } catch (Exception e) { + // TODO we should not catch Exception, and we should identify the offending additional equation + throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_ADDITIONAL_EQUATION); + } + } + + /** Return the dimension of one added set of equations. + * @param addEqu Additional equations used as a reference for selection + * @return dimension of the added set of equations + * @throws IllegalArgumentException if additional equation is unknown + */ + public int getAddedSetDimension(final AdditionalEquations addEqu) { + return selectStateAndEquations(addEqu).getAdditionalEquations().getDimension(); + } + + /** Get the current time derivative of the total state vector. + * @param t current value of the independent time variable + * @param y array containing the current value of the state vector + * @param yDot placeholder array where to put the time derivative of the state vector + */ + public void computeDerivatives(final double t, final double[] y, final double[] yDot) { + + // Add contribution for the main set of equations + int index = getMainSetDimension(); + double[] m = new double[index]; + double[] mDot = new double[index]; + // update current main state + System.arraycopy(y, 0, m, 0, index); + // compute derivatives + mainSet.computeDerivatives(t, m, mDot); + // update main state contribution in global array + System.arraycopy(mDot, 0, yDot, 0, index); + + // Add contribution for additional equations + for (final AdditionalStateAndEquations stateAndEqu : addedSets) { + final double[] p = stateAndEqu.getAdditionalState(); + final double[] pDot = stateAndEqu.getAdditionalStateDot(); + + // update current additional state + System.arraycopy(y, index, p, 0, p.length); + + // compute additional derivatives + stateAndEqu.getAdditionalEquations().computeDerivatives(t, m, mDot, p, pDot); + + // update each additional state contribution in global array + System.arraycopy(pDot, 0, yDot, index, p.length); + + // incrementing index + index += p.length; + } + + } + + /** Add a set of user-specified equations to be integrated along with + * the main set of equations. + * + * @param addEqu additional equations + * @see #setInitialAdditionalState(double[], AdditionalEquations) + * @see #getCurrentAdditionalState(AdditionalEquations) + */ + public void addAdditionalEquations(final AdditionalEquations addEqu) { + addedSets.add(new AdditionalStateAndEquations(addEqu)); + } + + /** Get the instance of the main set of equations. + * @return current value of the main set of equations. + */ + public FirstOrderDifferentialEquations getMainSet() { + return mainSet; + } + + /** Set initial additional state. + * @param addState additional state + * @param addEqu additional equations used as a reference for selection + * @throws IllegalArgumentException if additional equation is unknown + */ + public void setInitialAdditionalState(final double[] addState, final AdditionalEquations addEqu) { + selectStateAndEquations(addEqu).setAdditionalState(addState); + } + + /** Set current additional state. + *

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

+ * @param currentState total current state + * @throws IllegalArgumentException if additional equation is unknown + */ + public void setCurrentAdditionalState(final double[] currentState) + throws IllegalArgumentException { + int index = getMainSetDimension(); + try { + for (AdditionalStateAndEquations stateAndEqu : addedSets) { + final int addDim = stateAndEqu.getAdditionalEquations().getDimension(); + final double[] addState = new double[addDim]; + System.arraycopy(currentState, index, addState, 0, addDim); + stateAndEqu.setAdditionalState(addState); + index += addDim; + } + } catch (Exception e) { + // TODO we should not catch Exception, and we should identify the offending additional equation + throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_ADDITIONAL_EQUATION); + } + } + + /** Get current additional state. + * @param addEqu additional equations used as a reference for selection + * @return current additional state + * @throws IllegalArgumentException if additional equation is unknown + */ + public double[] getCurrentAdditionalState(final AdditionalEquations addEqu) { + return selectStateAndEquations(addEqu).getAdditionalState(); + } + + /** Get all current additional states accumulated. + * @return current additional states + * @throws IllegalArgumentException if additional equation is unknown + */ + public double[] getCurrentAdditionalStates() + throws IllegalArgumentException { + int index = 0; + final double[] cumulState = new double[getAddedSetsDimension()]; + try { + for (AdditionalStateAndEquations stateAndEqu : addedSets) { + final int addDim = stateAndEqu.getAdditionalEquations().getDimension(); + final double[] addState = stateAndEqu.getAdditionalState(); + System.arraycopy(addState, 0, cumulState, index, addDim); + index += addDim; + } + return cumulState; + } catch (Exception e) { + // TODO we should not catch Exception, and we should identify the offending additional equation + throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_ADDITIONAL_EQUATION); + } + } + + /** Select additional state and equations pair in the list. + * @param addEqu Additional equations used as a reference for selection + * @return additional state and equations pair + * @throws IllegalArgumentException if additional equation is unknown + */ + private AdditionalStateAndEquations selectStateAndEquations(final AdditionalEquations addEqu) + throws IllegalArgumentException { + for (AdditionalStateAndEquations stateAndEqu : addedSets) { + if (stateAndEqu.getAdditionalEquations() == addEqu) { + return stateAndEqu; + } + } + // TODO we should not catch Exception, and we should identify the offending additional equation + throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_ADDITIONAL_EQUATION); + } + +} Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderDifferentialEquations.java ------------------------------------------------------------------------------ svn:eol-style = native Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderDifferentialEquations.java ------------------------------------------------------------------------------ svn:keywords = "Author Date Id Revision" Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderIntegrator.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderIntegrator.java?rev=1175409&view=auto ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderIntegrator.java (added) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderIntegrator.java Sun Sep 25 15:04:39 2011 @@ -0,0 +1,65 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math.ode; + +import org.apache.commons.math.exception.MathIllegalArgumentException; +import org.apache.commons.math.exception.MathIllegalStateException; + +/** + * This interface represents a first order integrator for expandable + * differential equations. + *

+ * The classes devoted to solve expandable first order differential equations + * should implement this interface. The problems which can be handled should + * implement the {@link ExpandableFirstOrderDifferentialEquations} interface. + *

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

This method solves an Initial Value Problem (IVP).

+ *

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

+ *

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

+ * @param equations complete set of differential equations to integrate + * @param t0 initial time + * @param y0 initial value of the main state vector at t0 + * @param t target time for the integration + * (can be set to a value smaller than t0 for backward integration) + * @param y placeholder where to put the main state vector at each successful + * step (and hence at the end of integration), can be the same object as y0 + * @return stop time, will be the same as target time if integration reached its + * target, but may be different if some {@link + * org.apache.commons.math.ode.events.EventHandler} stops it at some point. + * @throws MathIllegalStateException if the integrator cannot perform integration + * @throws MathIllegalArgumentException if integration parameters are wrong (typically + * too small integration span) + */ + double integrate(ExpandableFirstOrderDifferentialEquations equations, + double t0, double[] y0, double t, double[] y) + throws MathIllegalStateException, MathIllegalArgumentException; + +} Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderIntegrator.java ------------------------------------------------------------------------------ svn:eol-style = native Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderIntegrator.java ------------------------------------------------------------------------------ svn:keywords = "Author Date Id Revision" Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/JacobianMatrices.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/JacobianMatrices.java?rev=1175409&view=auto ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/JacobianMatrices.java (added) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/JacobianMatrices.java Sun Sep 25 15:04:39 2011 @@ -0,0 +1,471 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math.ode; + +import java.lang.reflect.Array; +import java.util.ArrayList; +import java.util.List; + +import org.apache.commons.math.exception.DimensionMismatchException; +import org.apache.commons.math.exception.MathIllegalArgumentException; +import org.apache.commons.math.exception.util.LocalizedFormats; + +/** + * This class defines a set of {@link AdditionalEquations additional equations} to + * compute the jacobian matrices with respect to the initial state vector and, if + * any, to some parameters of the main ODE set. + *

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

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

+ * + * @see ExpandableFirstOrderDifferentialEquations + * @see FirstOrderDifferentialEquations + * @see MainStateJacobianProvider + * @see ParameterJacobianProvider + * @see ParameterizedODE + * + * @version $Id$ + * @since 3.0 + */ +public class JacobianMatrices implements AdditionalEquations { + + /** Expandable first order differential equation. */ + private ExpandableFirstOrderDifferentialEquations efode; + + /** FODE without exact main jacobian computation skill. */ + private FirstOrderDifferentialEquations fode = null; + + /** FODE with exact main jacobian computation skill. */ + private MainStateJacobianProvider jode = null; + + /** FODE without exact parameter jacobian computation skill. */ + private ParameterizedODE pode = null; + + /** FODE with exact parameter jacobian computation skill. */ + private List pjp = new ArrayList();; + + /** List of parameters selected for parameter jacobian computation. */ + private List selectedParameters = null; + + /** Main state vector dimension. */ + private int stateDim; + + /** Parameters dimension. */ + private int paramDim = 0; + + /** Current main state jacobian matrix in a row. */ + private double[] mainJacobianInARow; + + /** Current parameters jacobian matrices in a row. */ + private double[] parameterJacobiansInARow = null; + + /** Step used for finite difference computation of jacobian matrix + * w.r.t. main state vector. */ + private double[] hY = null; + + /** Boolean for fode consistency. */ + private boolean dirtyMainState = false; + + /** Boolean for selected parameters consistency. */ + private boolean dirtyParameter = false; + + /** Simple constructor for an additional equations set computing jacobian matrices. + *

This additional equations set is added internally to the expandable + * first order differential equations set thanks to the + * {@link ExpandableFirstOrderDifferentialEquations#addAdditionalEquations(AdditionalEquations)} + * method. + * @param extended the expandable first order differential equations set + * @param jode the main first order differential equations set to extend + * @exception IllegalArgumentException if jode does not match the main set to be extended given by + * {@link ExpandableFirstOrderDifferentialEquations#getMainSet() extended.getMainSet()} + */ + public JacobianMatrices(final ExpandableFirstOrderDifferentialEquations extended, + final MainStateJacobianProvider jode) + throws IllegalArgumentException { + + checkCompatibility(extended, jode); + + efode = extended; + stateDim = efode.getMainSetDimension(); + mainJacobianInARow = new double[stateDim * stateDim]; + this.jode = jode; + efode.addAdditionalEquations(this); + setInitialMainStateJacobian(); + } + + /** Simple constructor for an additional equations set computing jacobian matrices. + *

This additional equations set is added internally to the expandable + * first order differential equations set thanks to the + * {@link ExpandableFirstOrderDifferentialEquations#addAdditionalEquations(AdditionalEquations)} + * method. + * @param extended the expandable first order differential equations set + * @param fode the main first order differential equations set to extend + * @exception IllegalArgumentException if fode does not match the main set to be extended given by + * {@link ExpandableFirstOrderDifferentialEquations#getMainSet() extended.getMainSet()} + */ + public JacobianMatrices(final ExpandableFirstOrderDifferentialEquations extended, + final FirstOrderDifferentialEquations fode) + throws IllegalArgumentException { + + checkCompatibility(extended, fode); + + efode = extended; + stateDim = efode.getMainSetDimension(); + mainJacobianInARow = new double[stateDim * stateDim]; + this.fode = fode; + dirtyMainState = true; + efode.addAdditionalEquations(this); + setInitialMainStateJacobian(); + } + + /** Add a parameter jacobian provider. + * @param pjp the parameter jacobian provider to compute exactly the parameter jacobian matrix + */ + public void setParameterJacobianProvider(final ParameterJacobianProvider pjp) { + this.pjp.add(pjp); + } + + /** Add a parameter jacobian provider. + * @param pjp the parameterized ODE to compute by finite difference the parameter jacobian matrix + */ + public void setParameterizedODE(final ParameterizedODE pode) { + this.pode = pode; + dirtyParameter = true; + } + + /** Select the parameters to consider for jacobian matrices processing. + *

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

+ *

Note that each selection clears the previous selected parameters.

+ * + * @param parameters parameters to consider for jacobian matrices processing + * @exception IllegalArgumentException if one parameter is not supported + */ + public void selectParameters(final String... parameters) throws IllegalArgumentException { + + selectedParameters = new ArrayList(); + for (String param : parameters) { + selectedParameters.add(new ParameterConfiguration(param, Double.NaN)); + } + paramDim = parameters.length; + parameterJacobiansInARow = new double[paramDim * stateDim]; + setInitialParameterJacobians(); + } + + /** Set the step associated to a parameter in order to compute by finite + * difference the jacobian matrix. + *

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

+ *

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

+ *

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

+ * @param parameter parameter to consider for jacobian processing + * @param hP step for jacobian finite difference computation w.r.t. the specified parameter + * @see ParameterizedODE + * @exception IllegalArgumentException if the parameter is not supported + */ + public void setParameterStep(final String parameter, final double hP) { + + boolean found = false; + for (ParameterConfiguration param: selectedParameters) { + if (parameter.equals(param.getParameterName())) { + param.setHP(hP); + found = true; + dirtyParameter = true; + break; + } + } + if (!found) { + throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER, + parameter); + } + } + + /** Set the steps in order to compute by finite difference the jacobian + * matrix with respect to main state. + *

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

+ *

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

+ * @param hY step used for finite difference computation with respect to state vector + * @exception IllegalArgumentException if the hY has not the dimension of the main state + * given by {@link ExpandableFirstOrderDifferentialEquations#getMainSetDimension()} + */ + public void setMainStateSteps(final double[] hY) { + + if (fode != null) { + // Check dimension + checkDimension(stateDim, hY); + this.hY = hY.clone(); + dirtyMainState = true; + } + } + + /** Set the initial value of the jacobian matrix with respect to state. + * @param dYdY0 initial jacobian matrix w.r.t. state + * @exception IllegalArgumentException if matrix dimensions are incorrect + */ + public void setInitialMainStateJacobian(final double[][] dYdY0) { + + // Check dimensions + checkDimension(stateDim, dYdY0); + checkDimension(stateDim, dYdY0[0]); + + // store the matrix in row major order as a single dimension array + int index = 0; + for (final double[] row : dYdY0) { + System.arraycopy(row, 0, mainJacobianInARow, index, stateDim); + index += stateDim; + } + // set initial additional state value in expandable fode + efode.setInitialAdditionalState(mainJacobianInARow, this); + } + + /** Set the initial value of the jacobian matrix with respect to one parameter. + *

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

+ * @param pName parameter name + * @param dYdP initial jacobian matrix w.r.t. the parameter + * @exception IllegalArgumentException if matrix dimensions are incorrect + */ + public void setInitialParameterJacobian(final String pName, final double[] dYdP) { + + // Check dimensions + checkDimension(stateDim, dYdP); + + // store the matrix in a global single dimension array + boolean found = false; + int index = 0; + for (ParameterConfiguration param: selectedParameters) { + if (pName.equals(param.getParameterName())) { + System.arraycopy(dYdP, 0, parameterJacobiansInARow, index, stateDim); + double[] p = new double[this.getDimension()]; + index = stateDim * stateDim; + System.arraycopy(mainJacobianInARow, 0, p, 0, index); + System.arraycopy(parameterJacobiansInARow, 0, p, index, stateDim * paramDim); + // set initial additional state value in expandable fode + efode.setInitialAdditionalState(p, this); + found = true; + break; + } + index += stateDim; + } + if (! found) { + throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER, + pName); + } + } + + /** Set the default initial value of the jacobian matrix with respect to state. + *

dYdY0 is set to the identity matrix.

+ */ + public void setInitialMainStateJacobian() { + final double[][] dYdY0 = new double[stateDim][stateDim]; + for (int i = 0; i < stateDim; ++i) { + dYdY0[i][i] = 1.0; + } + setInitialMainStateJacobian(dYdY0); + } + + /** Set the default initial value of the jacobian matrix with respect to one parameter. + *

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

+ *

dYdP is set to the null matrix.

+ * @param pName parameter name + */ + public void setInitialParameterJacobian(final String pName) { + setInitialParameterJacobian(pName, new double[stateDim]); + } + + /** Set the default initial values of jacobian matrices with respect to all parameters. + */ + public void setInitialParameterJacobians() { + for (ParameterConfiguration param: selectedParameters) { + setInitialParameterJacobian(param.getParameterName()); + } + } + + /** Set default initial values for jacobian matrices. + *

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

+ */ + public void setInitialJacobians() { + setInitialMainStateJacobian(); + setInitialParameterJacobians(); + } + + /** Get the current value of the jacobian matrix with respect to state. + * @param dYdY0 current jacobian matrix with respect to state. + */ + public void getCurrentMainSetJacobian(final double[][] dYdY0) { + + // get current state for this set of equations from the expandable fode + double[] p = efode.getCurrentAdditionalState(this); + + int index = 0; + for (int i = 0; i < stateDim; i++) { + System.arraycopy(p, index, dYdY0[i], 0, stateDim); + index += stateDim; + } + + } + + /** Get the current value of the jacobian matrix with respect to one parameter. + * @param pName name of the parameter for the computed jacobian matrix + * @param dYdP current jacobian matrix with respect to the named parameter + */ + public void getCurrentParameterJacobian(String pName, final double[] dYdP) { + + // get current state for this set of equations from the expandable fode + double[] p = efode.getCurrentAdditionalState(this); + + int index = stateDim * stateDim; + for (ParameterConfiguration param: selectedParameters) { + if (param.getParameterName().equals(pName)) { + System.arraycopy(p, index, dYdP, 0, stateDim); + break; + } + index += stateDim; + } + + } + + /** {@inheritDoc} */ + public int getDimension() { + return stateDim * (stateDim + paramDim); + } + + /** {@inheritDoc} */ + public void computeDerivatives(final double t, final double[] y, final double[] yDot, + final double[] z, final double[] zDot) { + + // Lazy initialization + if (dirtyMainState) { + jode = new MainStateJacobianWrapper(fode, hY); + dirtyMainState = false; + } + if (dirtyParameter && (paramDim != 0)) { + pjp.add(new ParameterJacobianWrapper(jode, pode, selectedParameters)); + dirtyParameter = false; + } + + // variational equations: + // from d[dy/dt]/dy0 and d[dy/dt]/dp to d[dy/dy0]/dt and d[dy/dp]/dt + + // compute jacobian matrix with respect to main state + double[][] dFdY = new double[stateDim][stateDim]; + jode.computeMainStateJacobian(t, y, yDot, dFdY); + + // Dispatch jacobian matrix in the compound additional state vector + for (int i = 0; i < stateDim; ++i) { + final double[] dFdYi = dFdY[i]; + for (int j = 0; j < stateDim; ++j) { + double s = 0; + final int startIndex = j; + int zIndex = startIndex; + for (int l = 0; l < stateDim; ++l) { + s += dFdYi[l] * z[zIndex]; + zIndex += stateDim; + } + zDot[startIndex + i * stateDim] = s; + } + } + + if (paramDim != 0) { + // compute jacobian matrices with respect to parameters + double[] dFdP = new double[stateDim]; + int startIndex = stateDim * stateDim; + for (ParameterConfiguration param: selectedParameters) { + boolean found = false; + for (ParameterJacobianProvider provider: pjp) { + if (provider.isSupported(param.getParameterName())) { + try { + provider.computeParameterJacobian(t, y, yDot, param.getParameterName(), dFdP); + for (int i = 0; i < stateDim; ++i) { + final double[] dFdYi = dFdY[i]; + int zIndex = startIndex; + double s = dFdP[i]; + for (int l = 0; l < stateDim; ++l) { + s += dFdYi[l] * z[zIndex]; + zIndex++; + } + zDot[startIndex + i] = s; + } + startIndex += stateDim; + found = true; + break; + } catch (IllegalArgumentException iae) { + } + } + } + if (! found) { + throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER, + param); + } + } + } + + } + + /** Check compatibility between the main set in the expandable ode and an ordinary ode. + * @param expended expandable ode containing a main set + * @param ode single ode to check + * @throws MathIllegalArgumentException if single ode doesn't match the main ode set in the extended ode + */ + private void checkCompatibility(final ExpandableFirstOrderDifferentialEquations extended, + final FirstOrderDifferentialEquations ode) + throws MathIllegalArgumentException { + + if (!(ode == extended.getMainSet())) { + throw new MathIllegalArgumentException(LocalizedFormats.UNMATCHED_ODE_IN_EXTENDED_SET); + } + } + + /** Check array dimensions. + * @param expected expected dimension + * @param array (may be null if expected is 0) + * @throws DimensionMismatchException if the array dimension does not match the expected one + */ + private void checkDimension(final int expected, final Object array) + throws DimensionMismatchException { + int arrayDimension = (array == null) ? 0 : Array.getLength(array); + if (arrayDimension != expected) { + throw new DimensionMismatchException(arrayDimension, expected); + } + } + +} + Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/JacobianMatrices.java ------------------------------------------------------------------------------ svn:eol-style = native Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/JacobianMatrices.java ------------------------------------------------------------------------------ svn:keywords = "Author Date Id Revision" Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MainStateJacobianProvider.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MainStateJacobianProvider.java?rev=1175409&view=auto ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MainStateJacobianProvider.java (added) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MainStateJacobianProvider.java Sun Sep 25 15:04:39 2011 @@ -0,0 +1,36 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math.ode; + +/** Interface expanding {@link FirstOrderDifferentialEquations first order + * differential equations} in order to compute exactly the main state jacobian + * matrix for {@link JacobianMatrices partial derivatives equations}. + * + * @version $Id$ + * @since 3.0 + */ +public interface MainStateJacobianProvider extends FirstOrderDifferentialEquations { + + /** Compute the jacobian matrix of ODE with respect to main state. + * @param t current value of the independent time variable + * @param y array containing the current value of the main state vector + * @param yDot array containing the current value of the time derivative of the main state vector + * @param dFdY placeholder array where to put the jacobian matrix of the ODE w.r.t. the main state vector + */ + void computeMainStateJacobian(double t, double[] y, double[] yDot, double[][] dFdY); + +} Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MainStateJacobianProvider.java ------------------------------------------------------------------------------ svn:eol-style = native Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MainStateJacobianProvider.java ------------------------------------------------------------------------------ svn:keywords = "Author Date Id Revision" Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MainStateJacobianWrapper.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MainStateJacobianWrapper.java?rev=1175409&view=auto ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MainStateJacobianWrapper.java (added) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MainStateJacobianWrapper.java Sun Sep 25 15:04:39 2011 @@ -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; + } + } + +} Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MainStateJacobianWrapper.java ------------------------------------------------------------------------------ svn:eol-style = native Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MainStateJacobianWrapper.java ------------------------------------------------------------------------------ svn:keywords = "Author Date Id Revision" Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java?rev=1175409&r1=1175408&r2=1175409&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java Sun Sep 25 15:04:39 2011 @@ -403,7 +403,7 @@ public abstract class MultistepIntegrato } /** 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 MultistepIntegrato return dimension; } - /** {@inheritDoc} */ - public int getMainSetDimension() { - return mainSetDimension; - } } } Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterConfiguration.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterConfiguration.java?rev=1175409&view=auto ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterConfiguration.java (added) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterConfiguration.java Sun Sep 25 15:04:39 2011 @@ -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; + } + +} Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterConfiguration.java ------------------------------------------------------------------------------ svn:eol-style = native Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterConfiguration.java ------------------------------------------------------------------------------ svn:keywords = "Author Date Id Revision" Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianProvider.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianProvider.java?rev=1175409&view=auto ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianProvider.java (added) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianProvider.java Sun Sep 25 15:04:39 2011 @@ -0,0 +1,43 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math.ode; + +import org.apache.commons.math.exception.MathIllegalArgumentException; + +/** Interface to compute exactly jacobian matrix for some parameter + * when computing {@link JacobianMatrices partial derivatives equations}. + * + * @version $Id$ + * @since 3.0 + */ +public interface ParameterJacobianProvider extends Parameterizable { + + /** Compute the jacobian matrix of ODE with respect to one parameter. + *

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

+ * @param t current value of the independent time variable + * @param y array containing the current value of the main state vector + * @param yDot array containing the current value of the time derivative + * of the main state vector + * @param paramName name of the parameter to consider + * @param dFdP placeholder array where to put the jacobian matrix of the + * ODE with respect to the parameter + * @throws MathIllegalArgumentException if the parameter is not supported + */ + void computeParameterJacobian(double t, double[] y, double[] yDot, + String paramName, double[] dFdP) + throws MathIllegalArgumentException; +} Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianProvider.java ------------------------------------------------------------------------------ svn:eol-style = native Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianProvider.java ------------------------------------------------------------------------------ svn:keywords = "Author Date Id Revision" Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianWrapper.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianWrapper.java?rev=1175409&view=auto ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianWrapper.java (added) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianWrapper.java Sun Sep 25 15:04:39 2011 @@ -0,0 +1,91 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math.ode; + +import java.util.Collection; +import java.util.HashMap; +import java.util.Map; + +/** Wrapper class to compute jacobian matrices by finite differences for ODE + * which do not compute them by themselves. + * + * @version $Id$ + * @since 3.0 + */ +class ParameterJacobianWrapper implements ParameterJacobianProvider { + + /** Main ODE set. */ + private final FirstOrderDifferentialEquations fode; + + /** Raw ODE without jacobian computation skill to be wrapped into a ParameterJacobianProvider. */ + private final ParameterizedODE pode; + + /** Steps for finite difference computation of the jacobian df/dp w.r.t. parameters. */ + private final Map hParam; + + /** Wrap a {@link ParameterizedODE} into a {@link ParameterJacobianProvider}. + * @param fode main first order differential equations set + * @param pode additional problem, without parametre jacobian computation skill + * @param paramsAndSteps parameters and steps to compute the jacobians df/dp + * @see JacobianMatrices#setParameterStep(String, double) + */ + public ParameterJacobianWrapper(final FirstOrderDifferentialEquations fode, + final ParameterizedODE pode, + final Collection paramsAndSteps) { + this.fode = fode; + this.pode = pode; + this.hParam = new HashMap(); + + // set up parameters for jacobian computation + for (final ParameterConfiguration param : paramsAndSteps) { + final String name = param.getParameterName(); + if (pode.isSupported(name)) { + hParam.put(name, param.getHP()); + } + } + } + + /** {@inheritDoc} */ + public Collection getParametersNames() { + return pode.getParametersNames(); + } + + /** {@inheritDoc} */ + public boolean isSupported(String name) { + return pode.isSupported(name); + } + + /** {@inheritDoc} */ + public void computeParameterJacobian(double t, double[] y, double[] yDot, + String paramName, double[] dFdP) { + + final int n = fode.getDimension(); + final double[] tmpDot = new double[n]; + + // compute the jacobian df/dp w.r.t. parameter + final double p = pode.getParameter(paramName); + final double hP = hParam.get(paramName); + pode.setParameter(paramName, p + hP); + fode.computeDerivatives(t, y, tmpDot); + for (int i = 0; i < n; ++i) { + dFdP[i] = (tmpDot[i] - yDot[i]) / hP; + } + pode.setParameter(paramName, p); + + } + +} Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianWrapper.java ------------------------------------------------------------------------------ svn:eol-style = native Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianWrapper.java ------------------------------------------------------------------------------ svn:keywords = "Author Date Id Revision" Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/Parameterizable.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/Parameterizable.java?rev=1175409&view=auto ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/Parameterizable.java (added) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/Parameterizable.java Sun Sep 25 15:04:39 2011 @@ -0,0 +1,43 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math.ode; + +import java.util.Collection; + +/** This interface enables to process any parameterizable object. + * + * @version $Id$ + * @since 3.0 + */ + +public interface Parameterizable { + + /** Get the names of the supported parameters. + * @return parameters names + * @see #isSupported(String) + */ + Collection getParametersNames(); + + /** Check if a parameter is supported. + *

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

+ * @param name parameter name to check + * @return true if the parameter is supported + * @see #getParametersNames() + */ + boolean isSupported(String name); + +} Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/Parameterizable.java ------------------------------------------------------------------------------ svn:eol-style = native Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/Parameterizable.java ------------------------------------------------------------------------------ svn:keywords = "Author Date Id Revision" Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedODE.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedODE.java?rev=1175409&view=auto ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedODE.java (added) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedODE.java Sun Sep 25 15:04:39 2011 @@ -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; + +} Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedODE.java ------------------------------------------------------------------------------ svn:eol-style = native Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedODE.java ------------------------------------------------------------------------------ svn:keywords = "Author Date Id Revision" Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedWrapper.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedWrapper.java?rev=1175409&view=auto ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedWrapper.java (added) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedWrapper.java Sun Sep 25 15:04:39 2011 @@ -0,0 +1,77 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math.ode; + +import java.util.ArrayList; +import java.util.Collection; + +import org.apache.commons.math.exception.MathIllegalArgumentException; +import org.apache.commons.math.exception.util.LocalizedFormats; + + +/** Wrapper class enabling {@link FirstOrderDifferentialEquations basic simple} + * ODE instances to be used when processing {@link JacobianMatrices}. + * + * @version $Id$ + * @since 3.0 + */ +class ParameterizedWrapper implements ParameterizedODE { + + /** Basic FODE without parameter. */ + private final FirstOrderDifferentialEquations fode; + + /** Simple constructor. + * @param ode original first order differential equations + */ + public ParameterizedWrapper(final FirstOrderDifferentialEquations ode) { + this.fode = ode; + } + + /** {@inheritDoc} */ + public int getDimension() { + return fode.getDimension(); + } + + /** {@inheritDoc} */ + public void computeDerivatives(double t, double[] y, double[] yDot) { + fode.computeDerivatives(t, y, yDot); + } + + /** {@inheritDoc} */ + public Collection getParametersNames() { + return new ArrayList(); + } + + /** {@inheritDoc} */ + public boolean isSupported(String name) { + return false; + } + + /** {@inheritDoc} */ + public double getParameter(String name) + throws MathIllegalArgumentException { + if (!isSupported(name)) { + throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER, name); + } + return Double.NaN; + } + + /** {@inheritDoc} */ + public void setParameter(String name, double value) { + } + +} Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedWrapper.java ------------------------------------------------------------------------------ svn:eol-style = native Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedWrapper.java ------------------------------------------------------------------------------ svn:keywords = "Author Date Id Revision" Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java?rev=1175409&r1=1175408&r2=1175409&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java Sun Sep 25 15:04:39 2011 @@ -20,7 +20,7 @@ package org.apache.commons.math.ode.nons 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 ex /** {@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 ex } 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;