You are viewing a plain text version of this content. The canonical link for it is here.
Posted to dev@commons.apache.org by lu...@free.fr on 2010/03/01 18:23:59 UTC
[math] adding a way to compute derivatives wrt initial states and
parameters in ODE (was Re: svn commit: r917592)
Hi,
This commit implements a new feature I wanted to have for a long time. It sits on top of ODE integrators which solve initial value problems of the form:
y(t0) = y0
y'(t) = f(t, y)
In addition to compute y(t) at any time, it also computes the derivatives (jacobians) dy/dy0 and dy/dp where y0 and p are respectively the initial state and some problem parameters. This is useful for example when you need to adjust the initial state to match some conditions, like in Boundary Value Problems or optimization.
This is a first (working) implementation and it needs some polishing. I would like to have some feedback about it as I work on it. My first dumb question is: could someone find better names for the classes and interfaces I created ? As you know, I use rather long names but ParameterizedFirstOrderDifferentialEquationsWithPartials is probably too much ...
My second question is: how should I implement all the bells and whistles the ODE pacakge has, essentially for discrete events handling and step handling ? I was not able to set up the integrator using inheritance, should I try harder to force this problem into the existing ODE hierarchy or should I keep duplicate the hierarchy and add the parameters for jacobians (the existing hierarchy has no way to pass jacobians information in the methods signatures) ?
This feature is a very important one for many applications and it will open the way to a complete solver for boundary value problems, with simple or multiple shooting methods.
Thanks in advance for your advices
Luc
----- luc@apache.org a écrit :
> Author: luc
> Date: Mon Mar 1 17:07:47 2010
> New Revision: 917592
>
> URL: http://svn.apache.org/viewvc?rev=917592&view=rev
> Log:
> Added a way to compute both the final state in an Initial Value
> Problem (IVP)
> for Ordinary Differential Equations (ODE) and its derivatives with
> respect to
> initial state and with respect to some problem parameters. This allows
> wrapping
> ODE solvers into optimization or root finding algorithms, which in
> turn can be
> used to solve Boundary Value Problems (BVP). There are no
> implementations (yet)
> of BVP solvers in the library.
>
> Added:
>
> commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/EnhancedFirstOrderIntegrator.java
> (with props)
>
> commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedFirstOrderDifferentialEquations.java
> (with props)
>
> commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedFirstOrderDifferentialEquationsWithPartials.java
> (with props)
>
> commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/EnhancedFirstOrderIntegratorTest.java
> (with props)
> Modified:
> commons/proper/math/trunk/src/site/xdoc/changes.xml
>
> Added:
> commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/EnhancedFirstOrderIntegrator.java
> URL:
> http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/EnhancedFirstOrderIntegrator.java?rev=917592&view=auto
> ==============================================================================
> ---
> commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/EnhancedFirstOrderIntegrator.java
> (added)
> +++
> commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/EnhancedFirstOrderIntegrator.java
> Mon Mar 1 17:07:47 2010
> @@ -0,0 +1,325 @@
> +/*
> + * 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 org.apache.commons.math.MathRuntimeException;
> +
> +/** This class enhances a first order integrator for differential
> equations to
> + * compute also partial derivatives of the solution with respect to
> initial state
> + * and parameters.
> + * <p>In order to compute both the state and its derivatives, the ODE
> problem
> + * is extended with jacobians of the raw ODE and the variational
> equations are
> + * added to form a new compound problem of higher dimension. If the
> original ODE
> + * problem has dimension n and there are p parameters, the compound
> problem will
> + * have dimension n × (1 + n + k).</p>
> + * @see ParameterizedFirstOrderDifferentialEquations
> + * @see ParameterizedFirstOrderDifferentialEquationsWithPartials
> + * @version $Revision$ $Date$
> + * @since 2.1
> + */
> +public class EnhancedFirstOrderIntegrator {
> +
> + /** Underlying integrator for compound problem. */
> + private final FirstOrderIntegrator integrator;
> +
> + /** Raw equations to integrate. */
> + private final
> ParameterizedFirstOrderDifferentialEquationsWithPartials ode;
> +
> + /** Build an enhanced integrator using internal differentiation
> to compute jacobians.
> + * @param integrator underlying integrator to solve the compound
> problem
> + * @param ode original problem (f in the equation y' = f(t, y))
> + * @param p parameters array (may be null if {@link
> + *
> ParameterizedFirstOrderDifferentialEquations#getParametersDimension()
> + * getParametersDimension()} from original problem is zero)
> + * @param hY step sizes to use for computing the jacobian df/dy,
> must have the
> + * same dimension as the original problem
> + * @param hP step sizes to use for computing the jacobian df/dp,
> must have the
> + * same dimension as the original problem parameters dimension
> + * @see #EnhancedFirstOrderIntegrator(FirstOrderIntegrator,
> + * ParameterizedFirstOrderDifferentialEquationsWithPartials)
> + */
> + public EnhancedFirstOrderIntegrator(final FirstOrderIntegrator
> integrator,
> + final
> ParameterizedFirstOrderDifferentialEquations ode,
> + final double[] p, final
> double[] hY, final double[] hP) {
> + checkDimension(ode.getDimension(), hY);
> + checkDimension(ode.getParametersDimension(), p);
> + checkDimension(ode.getParametersDimension(), hP);
> + this.integrator = integrator;
> + this.ode = new FiniteDifferencesWrapper(ode, p, hY, hP);
> + }
> +
> + /** Build an enhanced integrator using ODE builtin jacobian
> computation features.
> + * @param integrator underlying integrator to solve the compound
> problem
> + * @param ode original problem, which can compute the jacobians
> by itself
> + * @see #EnhancedFirstOrderIntegrator(FirstOrderIntegrator,
> + * ParameterizedFirstOrderDifferentialEquations, double[],
> double[], double[])
> + */
> + public EnhancedFirstOrderIntegrator(final FirstOrderIntegrator
> integrator,
> + final
> ParameterizedFirstOrderDifferentialEquationsWithPartials ode) {
> + this.integrator = integrator;
> + this.ode = ode;
> + }
> +
> + /** Integrate the differential equations and the variational
> equations up to the given time.
> + * <p>This method solves an Initial Value Problem (IVP) and also
> computes the derivatives
> + * of the solution with respect to initial state and parameters.
> This can be used as
> + * a basis to solve Boundary Value Problems (BVP).</p>
> + * <p>Since this method stores some internal state variables
> made
> + * available in its public interface during integration ({@link
> + * #getCurrentSignedStepsize()}), it is <em>not</em>
> thread-safe.</p>
> + * @param equations differential equations to integrate
> + * @param t0 initial time
> + * @param y0 initial value of the state vector at t0
> + * @param dY0dP initial value of the state vector derivative with
> respect to the
> + * parameters at t0
> + * @param t target time for the integration
> + * (can be set to a value smaller than <code>t0</code> for
> backward integration)
> + * @param y placeholder where to put the state vector at each
> successful
> + * step (and hence at the end of integration), can be the same
> object as y0
> + * @param dYdY0 placeholder where to put the state vector
> derivative with respect
> + * to the initial state (dy[i]/dy0[j] is in element array
> dYdY0[i][j]) at each successful
> + * step (and hence at the end of integration)
> + * @param dYdP placeholder where to put the state vector
> derivative with respect
> + * to the parameters (dy[i]/dp[j] is in element array dYdP[i][j])
> at each successful
> + * step (and hence at the end of integration)
> + * @return stop time, will be the same as target time if
> integration reached its
> + * target, but may be different if some event handler stops it at
> some point.
> + * @throws IntegratorException if the integrator cannot perform
> integration
> + * @throws DerivativeException this exception is propagated to
> the caller if
> + * the underlying user function triggers one
> + */
> + public double integrate(final double t0, final double[] y0, final
> double[][] dY0dP,
> + final double t, final double[] y,
> + final double[][] dYdY0, final double[][]
> dYdP)
> + throws DerivativeException, IntegratorException {
> +
> + final int n = ode.getDimension();
> + final int k = ode.getParametersDimension();
> + checkDimension(n, y0);
> + checkDimension(n, y);
> + checkDimension(n, dYdY0);
> + checkDimension(n, dYdY0[0]);
> + if (k != 0) {
> + checkDimension(n, dY0dP);
> + checkDimension(k, dY0dP[0]);
> + checkDimension(n, dYdP);
> + checkDimension(k, dYdP[0]);
> + }
> +
> + // the compound state z contains the raw state y and its
> derivatives
> + // with respect to initial state y0 and to parameters p
> + // y[i] is stored in z[i]
> + // dy[i]/dy0[j] is stored in z[n + i * n + j]
> + // dy[i]/dp[j] is stored in z[n * (n + 1) + i * k + j]
> + final int q = n * (1 + n + k);
> +
> + // set up initial state, including partial derivatives
> + final double[] z = new double[q];
> + System.arraycopy(y0, 0, z, 0, n);
> + for (int i = 0; i < n; ++i) {
> +
> + // set diagonal element of dy/dy0 to 1.0 at t = t0
> + z[i * (1 + n) + n] = 1.0;
> +
> + // set initial derivatives with respect to parameters
> + System.arraycopy(dY0dP[i], 0, z, n * (n + 1) + i * k,
> k);
> +
> + }
> +
> + // integrate the compound state variational equations
> + final double stopTime = integrator.integrate(new
> FirstOrderDifferentialEquations() {
> +
> + /** Current state. */
> + private final double[] y = new double[n];
> +
> + /** Time derivative of the current state. */
> + private final double[] yDot = new double[n];
> +
> + /** Derivatives of yDot with respect to state. */
> + private final double[][] dFdY = new double[n][n];
> +
> + /** Derivatives of yDot with respect to parameters. */
> + private final double[][] dFdP = new double[n][k];
> +
> + /** {@inheritDoc} */
> + public int getDimension() {
> + return q;
> + }
> +
> + /** {@inheritDoc} */
> + public void computeDerivatives(final double t, final
> double[] z, final double[] zDot)
> + throws DerivativeException {
> +
> + // compute raw ODE and its jacobians: dy/dt,
> d[dy/dt]/dy0 and d[dy/dt]/dp
> + System.arraycopy(z, 0, y, 0, n);
> + ode.computeDerivatives(t, y, yDot);
> + ode.computeJacobians(t, y, yDot, dFdY, dFdP);
> +
> + // state part of the compound equations
> + System.arraycopy(yDot, 0, zDot, 0, n);
> +
> + // variational equations: from d[dy/dt]/dy0 to
> d[dy/dy0]/dt
> + for (int i = 0; i < n; ++i) {
> + final double[] dFdYi = dFdY[i];
> + for (int j = 0; j < n; ++j) {
> + double s = 0;
> + int zIndex = n + j;
> + for (int l = 0; l < n; ++l) {
> + s += dFdYi[l] * z[zIndex];
> + zIndex += l;
> + }
> + zDot[n + i * n + j] = s;
> + }
> + }
> +
> + // variational equations: d[dy/dt]/dy0 and
> d[dy/dt]/dp to d[dy/dp]/dt
> + for (int i = 0; i < n; ++i) {
> + final double[] dFdYi = dFdY[i];
> + final double[] dFdPi = dFdP[i];
> + for (int j = 0; j < k; ++j) {
> + double s = dFdPi[j];
> + int zIndex = n * (n + 1)+ j;
> + for (int l = 0; l < n; ++l) {
> + s += dFdYi[l] * z[zIndex];
> + zIndex += k;
> + }
> + zDot[n * (n + 1) + i * k + j] = s;
> + }
> + }
> +
> + }
> +
> + }, t0, z, t, z);
> +
> + // dispatch the final compound state into the state and
> partial derivatives arrays
> + System.arraycopy(z, 0, y, 0, n);
> + for (int i = 0; i < n; ++i) {
> + System.arraycopy(z, n * (i + 1), dYdY0[i], 0, n);
> + }
> + for (int i = 0; i < n; ++i) {
> + System.arraycopy(z, n * (n + 1) + i * k, dYdP[i], 0, k);
> + }
> +
> + return stopTime;
> +
> + }
> +
> + /** Check array dimensions.
> + * @param expected expected dimension
> + * @param array (may be null if expected is 0)
> + * @throws IllegalArgumentException if the array dimension does
> not match the expected one
> + */
> + private void checkDimension(final int expected, final Object
> array)
> + throws IllegalArgumentException {
> + int arrayDimension = (array == null) ? 0 :
> Array.getLength(array);
> + if (arrayDimension != expected) {
> + throw
> MathRuntimeException.createIllegalArgumentException(
> + "dimension mismatch {0} != {1}", arrayDimension,
> expected);
> + }
> + }
> +
> + /** Wrapper class to compute jacobians by finite differences for
> ODE which do not compute them themselves. */
> + private static class FiniteDifferencesWrapper
> + implements
> ParameterizedFirstOrderDifferentialEquationsWithPartials {
> +
> + /** Raw ODE without jacobians computation. */
> + private final ParameterizedFirstOrderDifferentialEquations
> ode;
> +
> + /** Parameters array (may be null if parameters dimension
> from original problem is zero) */
> + private final double[] p;
> +
> + /** Step sizes to use for computing the jacobian df/dy. */
> + private final double[] hY;
> +
> + /** Step sizes to use for computing the jacobian df/dp. */
> + private final double[] hP;
> +
> + /** Temporary array for state derivatives used to compute
> jacobians. */
> + private final double[] tmpDot;
> +
> + /** Simple constructor.
> + * @param ode original ODE problem, without jacobians
> computations
> + * @param p parameters array (may be null if parameters
> dimension from original problem is zero)
> + * @param hY step sizes to use for computing the jacobian
> df/dy
> + * @param hP step sizes to use for computing the jacobian
> df/dp
> + */
> + public FiniteDifferencesWrapper(final
> ParameterizedFirstOrderDifferentialEquations ode,
> + final double[] p, final
> double[] hY, final double[] hP) {
> + this.ode = ode;
> + this.p = p.clone();
> + this.hY = hY.clone();
> + this.hP = hP.clone();
> + tmpDot = new double[ode.getDimension()];
> + }
> +
> + /** {@inheritDoc} */
> + public int getDimension() {
> + return ode.getDimension();
> + }
> +
> + /** {@inheritDoc} */
> + public void computeDerivatives(double t, double[] y, double[]
> yDot) throws DerivativeException {
> + ode.computeDerivatives(t, y, yDot);
> + }
> +
> + /** {@inheritDoc} */
> + public int getParametersDimension() {
> + return ode.getParametersDimension();
> + }
> +
> + /** {@inheritDoc} */
> + public void setParameter(int i, double value) {
> + ode.setParameter(i, value);
> + }
> +
> + /** {@inheritDoc} */
> + public void computeJacobians(double t, double[] y, double[]
> yDot,
> + double[][] dFdY, double[][]
> dFdP)
> + throws DerivativeException {
> +
> + final int n = ode.getDimension();
> + final int k = ode.getParametersDimension();
> +
> + // compute df/dy where f is the ODE and y is the state
> array
> + for (int j = 0; j < n; ++j) {
> + final double savedYj = y[j];
> + y[j] += hY[j];
> + ode.computeDerivatives(t, y, tmpDot);
> + for (int i = 0; i < n; ++i) {
> + dFdY[i][j] = (tmpDot[i] - yDot[i]) / hY[j];
> + }
> + y[j] = savedYj;
> + }
> +
> + // compute df/dp where f is the ODE and p is the
> parameters array
> + for (int j = 0; j < k; ++j) {
> + ode.setParameter(j, p[j] + hP[j]);
> + ode.computeDerivatives(t, y, tmpDot);
> + for (int i = 0; i < n; ++i) {
> + dFdP[i][j] = (tmpDot[i] - yDot[i]) / hP[j];
> + }
> + ode.setParameter(j, p[j]);
> + }
> +
> + }
> +
> + }
> +
> +}
>
> Propchange:
> commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/EnhancedFirstOrderIntegrator.java
> ------------------------------------------------------------------------------
> svn:eol-style = native
>
> Propchange:
> commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/EnhancedFirstOrderIntegrator.java
> ------------------------------------------------------------------------------
> svn:keywords = Author Date Id Revision
>
> Added:
> commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedFirstOrderDifferentialEquations.java
> URL:
> http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedFirstOrderDifferentialEquations.java?rev=917592&view=auto
> ==============================================================================
> ---
> commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedFirstOrderDifferentialEquations.java
> (added)
> +++
> commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedFirstOrderDifferentialEquations.java
> Mon Mar 1 17:07:47 2010
> @@ -0,0 +1,45 @@
> +/*
> + * Licensed to the Apache Software Foundation (ASF) under one or
> more
> + * contributor license agreements. See the NOTICE file distributed
> with
> + * this work for additional information regarding copyright
> ownership.
> + * The ASF licenses this file to You under the Apache License,
> Version 2.0
> + * (the "License"); you may not use this file except in compliance
> with
> + * the License. You may obtain a copy of the License at
> + *
> + * http://www.apache.org/licenses/LICENSE-2.0
> + *
> + * Unless required by applicable law or agreed to in writing,
> software
> + * distributed under the License is distributed on an "AS IS" BASIS,
> + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
> implied.
> + * See the License for the specific language governing permissions
> and
> + * limitations under the License.
> + */
> +
> +package org.apache.commons.math.ode;
> +
> +
> +/** This interface represents {@link FirstOrderDifferentialEquations
> + * first order differential equations} with parameters.
> + *
> + * @see EnhancedFirstOrderIntegrator
> + *
> + * @version $Revision$ $Date$
> + * @since 2.1
> + */
> +
> +public interface ParameterizedFirstOrderDifferentialEquations
> + extends FirstOrderDifferentialEquations {
> +
> + /** Get the number of parameters.
> + * @return number of parameters
> + */
> + int getParametersDimension();
> +
> + /** Set a parameter.
> + * @param i index of the parameters (must be between 0
> + * and {@link #getParametersDimension() getParametersDimension()
> - 1})
> + * @param value value for the parameter
> + */
> + void setParameter(int i, double value);
> +
> +}
>
> Propchange:
> commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedFirstOrderDifferentialEquations.java
> ------------------------------------------------------------------------------
> svn:eol-style = native
>
> Propchange:
> commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedFirstOrderDifferentialEquations.java
> ------------------------------------------------------------------------------
> svn:keywords = Author Date Id Revision
>
> Added:
> commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedFirstOrderDifferentialEquationsWithPartials.java
> URL:
> http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedFirstOrderDifferentialEquationsWithPartials.java?rev=917592&view=auto
> ==============================================================================
> ---
> commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedFirstOrderDifferentialEquationsWithPartials.java
> (added)
> +++
> commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedFirstOrderDifferentialEquationsWithPartials.java
> Mon Mar 1 17:07:47 2010
> @@ -0,0 +1,45 @@
> +/*
> + * Licensed to the Apache Software Foundation (ASF) under one or
> more
> + * contributor license agreements. See the NOTICE file distributed
> with
> + * this work for additional information regarding copyright
> ownership.
> + * The ASF licenses this file to You under the Apache License,
> Version 2.0
> + * (the "License"); you may not use this file except in compliance
> with
> + * the License. You may obtain a copy of the License at
> + *
> + * http://www.apache.org/licenses/LICENSE-2.0
> + *
> + * Unless required by applicable law or agreed to in writing,
> software
> + * distributed under the License is distributed on an "AS IS" BASIS,
> + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
> implied.
> + * See the License for the specific language governing permissions
> and
> + * limitations under the License.
> + */
> +
> +package org.apache.commons.math.ode;
> +
> +
> +/** This interface represents {@link
> ParameterizedFirstOrderDifferentialEquations
> + * first order differential equations} with parameters and partial
> derivatives.
> + *
> + * @see EnhancedFirstOrderIntegrator
> + *
> + * @version $Revision$ $Date$
> + * @since 2.1
> + */
> +
> +public interface
> ParameterizedFirstOrderDifferentialEquationsWithPartials
> + extends ParameterizedFirstOrderDifferentialEquations {
> +
> + /** Compute the partial derivatives of ODE with respect to
> state.
> + * @param t current value of the independent <I>time</I>
> variable
> + * @param y array containing the current value of the state
> vector
> + * @param yDot array containing the current value of the time
> derivative of the state vector
> + * @param dFdY placeholder array where to put the jacobian of the
> ODE with respect to the state vector
> + * @param dFdP placeholder array where to put the jacobian of the
> ODE with respect to the parameters
> + * @throws DerivativeException this exception is propagated to
> the caller if the
> + * underlying user function triggers one
> + */
> + void computeJacobians(double t, double[] y, double[] yDot,
> double[][] dFdY, double[][] dFdP)
> + throws DerivativeException;
> +
> +}
>
> Propchange:
> commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedFirstOrderDifferentialEquationsWithPartials.java
> ------------------------------------------------------------------------------
> svn:eol-style = native
>
> Propchange:
> commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedFirstOrderDifferentialEquationsWithPartials.java
> ------------------------------------------------------------------------------
> svn:keywords = Author Date Id Revision
>
> Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
> URL:
> http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=917592&r1=917591&r2=917592&view=diff
> ==============================================================================
> --- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
> +++ commons/proper/math/trunk/src/site/xdoc/changes.xml Mon Mar 1
> 17:07:47 2010
> @@ -39,6 +39,14 @@
> </properties>
> <body>
> <release version="2.1" date="TBD" description="TBD">
> + <action dev="luc" type="add" >
> + Added a way to compute both the final state in an Initial
> Value Problem (IVP)
> + for Ordinary Differential Equations (ODE) and its derivatives
> with respect to
> + initial state and with respect to some problem parameters.
> This allows wrapping
> + ODE solvers into optimization or root finding algorithms,
> which in turn can be
> + used to solve Boundary Value Problems (BVP). There are no
> implementations (yet)
> + of BVP solvers in the library.
> + </action>
> <action dev="luc" type="fix" issue="MATH-344" >
> Fixed wrong return values when enpoints are roots in Brent
> solver with
> a user provided initial guess
>
> Added:
> commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/EnhancedFirstOrderIntegratorTest.java
> URL:
> http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/EnhancedFirstOrderIntegratorTest.java?rev=917592&view=auto
> ==============================================================================
> ---
> commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/EnhancedFirstOrderIntegratorTest.java
> (added)
> +++
> commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/EnhancedFirstOrderIntegratorTest.java
> Mon Mar 1 17:07:47 2010
> @@ -0,0 +1,179 @@
> +/*
> + * 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.junit.Assert;
> +import org.junit.Test;
> +
> +public class EnhancedFirstOrderIntegratorTest {
> +
> + @Test
> + public void testLowAccuracyExternalDifferentiation()
> + throws IntegratorException, DerivativeException {
> + FirstOrderIntegrator integ =
> + new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-4,
> 1.0e-4);
> + double hP = 1.0e-12;
> + SummaryStatistics residuals0 = new SummaryStatistics();
> + SummaryStatistics residuals1 = new SummaryStatistics();
> + for (double b = 2.88; b < 3.08; b += 0.001) {
> + Brusselator brusselator = new Brusselator(b);
> + double[] y = { 1.3, b };
> + integ.integrate(brusselator, 0, y, 20.0, y);
> + double[] yP = { 1.3, b + hP };
> + brusselator.setParameter(0, b + hP);
> + integ.integrate(brusselator, 0, yP, 20.0, yP);
> + residuals0.addValue((yP[0] - y[0]) / hP -
> brusselator.dYdP0());
> + residuals1.addValue((yP[1] - y[1]) / hP -
> brusselator.dYdP1());
> + }
> + Assert.assertTrue((residuals0.getMax() - residuals0.getMin())
> > 600);
> + Assert.assertTrue(residuals0.getStandardDeviation() > 30);
> + Assert.assertTrue((residuals1.getMax() - residuals1.getMin())
> > 800);
> + Assert.assertTrue(residuals1.getStandardDeviation() > 50);
> + }
> +
> + @Test
> + public void testHighAccuracyExternalDifferentiation()
> + throws IntegratorException, DerivativeException {
> + FirstOrderIntegrator integ =
> + new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-10,
> 1.0e-10);
> + double hP = 1.0e-12;
> + SummaryStatistics residuals0 = new SummaryStatistics();
> + SummaryStatistics residuals1 = new SummaryStatistics();
> + for (double b = 2.88; b < 3.08; b += 0.001) {
> + Brusselator brusselator = new Brusselator(b);
> + double[] y = { 1.3, b };
> + integ.integrate(brusselator, 0, y, 20.0, y);
> + double[] yP = { 1.3, b + hP };
> + brusselator.setParameter(0, b + hP);
> + integ.integrate(brusselator, 0, yP, 20.0, yP);
> + residuals0.addValue((yP[0] - y[0]) / hP -
> brusselator.dYdP0());
> + residuals1.addValue((yP[1] - y[1]) / hP -
> brusselator.dYdP1());
> + }
> + Assert.assertTrue((residuals0.getMax() - residuals0.getMin())
> > 0.02);
> + Assert.assertTrue((residuals0.getMax() - residuals0.getMin())
> < 0.03);
> + Assert.assertTrue(residuals0.getStandardDeviation() >
> 0.003);
> + Assert.assertTrue(residuals0.getStandardDeviation() <
> 0.004);
> + Assert.assertTrue((residuals1.getMax() - residuals1.getMin())
> > 0.04);
> + Assert.assertTrue((residuals1.getMax() - residuals1.getMin())
> < 0.05);
> + Assert.assertTrue(residuals1.getStandardDeviation() >
> 0.006);
> + Assert.assertTrue(residuals1.getStandardDeviation() <
> 0.007);
> + }
> +
> + @Test
> + public void testInternalDifferentiation()
> + throws IntegratorException, DerivativeException {
> + FirstOrderIntegrator integ =
> + new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-4,
> 1.0e-4);
> + double hP = 1.0e-12;
> + SummaryStatistics residuals0 = new SummaryStatistics();
> + SummaryStatistics residuals1 = new SummaryStatistics();
> + for (double b = 2.88; b < 3.08; b += 0.001) {
> + Brusselator brusselator = new Brusselator(b);
> + brusselator.setParameter(0, b);
> + double[] z = { 1.3, b };
> + double[][] dZdZ0 = new double[2][2];
> + double[][] dZdP = new double[2][1];
> + double hY = 1.0e-12;
> + EnhancedFirstOrderIntegrator extInt =
> + new EnhancedFirstOrderIntegrator(integ, brusselator,
> new double[] { b },
> + new double[] { hY,
> hY }, new double[] { hP });
> + extInt.integrate(0, z, new double[][] { { 0.0 }, { 1.0 }
> }, 20.0, z, dZdZ0, dZdP);
> + residuals0.addValue(dZdP[0][0] - brusselator.dYdP0());
> + residuals1.addValue(dZdP[1][0] - brusselator.dYdP1());
> + }
> + Assert.assertTrue((residuals0.getMax() - residuals0.getMin())
> < 0.006);
> + Assert.assertTrue(residuals0.getStandardDeviation() <
> 0.0009);
> + Assert.assertTrue((residuals1.getMax() - residuals1.getMin())
> < 0.006);
> + Assert.assertTrue(residuals1.getStandardDeviation() <
> 0.0012);
> + }
> +
> + @Test
> + public void testAnalyticalDifferentiation()
> + throws IntegratorException, DerivativeException {
> + FirstOrderIntegrator integ =
> + new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-4,
> 1.0e-4);
> + SummaryStatistics residuals0 = new SummaryStatistics();
> + SummaryStatistics residuals1 = new SummaryStatistics();
> + for (double b = 2.88; b < 3.08; b += 0.001) {
> + Brusselator brusselator = new Brusselator(b);
> + brusselator.setParameter(0, b);
> + double[] z = { 1.3, b };
> + double[][] dZdZ0 = new double[2][2];
> + double[][] dZdP = new double[2][1];
> + EnhancedFirstOrderIntegrator extInt =
> + new EnhancedFirstOrderIntegrator(integ,
> brusselator);
> + extInt.integrate(0, z, new double[][] { { 0.0 }, { 1.0 }
> }, 20.0, z, dZdZ0, dZdP);
> + residuals0.addValue(dZdP[0][0] - brusselator.dYdP0());
> + residuals1.addValue(dZdP[1][0] - brusselator.dYdP1());
> + }
> + Assert.assertTrue((residuals0.getMax() - residuals0.getMin())
> < 0.004);
> + Assert.assertTrue(residuals0.getStandardDeviation() <
> 0.001);
> + Assert.assertTrue((residuals1.getMax() - residuals1.getMin())
> < 0.005);
> + Assert.assertTrue(residuals1.getStandardDeviation() <
> 0.001);
> + }
> +
> + private static class Brusselator implements
> ParameterizedFirstOrderDifferentialEquationsWithPartials {
> +
> + private double b;
> +
> + public Brusselator(double b) {
> + this.b = b;
> + }
> +
> + public int getDimension() {
> + return 2;
> + }
> +
> + public void setParameter(int i, double p) {
> + b = p;
> + }
> +
> + public int getParametersDimension() {
> + return 1;
> + }
> +
> + public void computeDerivatives(double t, double[] y, double[]
> yDot) {
> + double prod = y[0] * y[0] * y[1];
> + yDot[0] = 1 + prod - (b + 1) * y[0];
> + yDot[1] = b * y[0] - prod;
> + }
> +
> + public void computeJacobians(double t, double[] y, double[]
> yDot, double[][] dFdY, double[][] dFdP) {
> + double p = 2 * y[0] * y[1];
> + double y02 = y[0] * y[0];
> + dFdY[0][0] = p - (1 + b);
> + dFdY[0][1] = y02;
> + dFdY[1][0] = b - p;
> + dFdY[1][1] = -y02;
> + dFdP[0][0] = -y[0];
> + dFdP[1][0] = y[0];
> + }
> +
> + public double dYdP0() {
> + return -1087.8787631970476 + (1050.4387741821572 +
> (-338.90621620263096 + 36.51793006801084 * b) * b) * b;
> + }
> +
> + public double dYdP1() {
> + return 1499.0904666097015 + (-1434.9574631810726 +
> (459.71079478756945 - 49.29949940968588 * b) * b) * b;
> + }
> +
> + };
> +
> +}
>
> Propchange:
> commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/EnhancedFirstOrderIntegratorTest.java
> ------------------------------------------------------------------------------
> svn:eol-style = native
>
> Propchange:
> commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/EnhancedFirstOrderIntegratorTest.java
> ------------------------------------------------------------------------------
> svn:keywords = Author Date Id Revision
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org