Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java?rev=1175409&r1=1175408&r2=1175409&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.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.MultistepIntegrator; @@ -86,7 +86,7 @@ public abstract class AdamsIntegrator ex /** {@inheritDoc} */ @Override - public abstract double integrate(final FirstOrderDifferentialEquations equations, + public abstract double integrate(final ExpandableFirstOrderDifferentialEquations equations, final double t0, final double[] y0, final double t, final double[] y) throws MathIllegalStateException, MathIllegalArgumentException;
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java?rev=1175409&r1=1175408&r2=1175409&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java Sun Sep 25 15:04:39 2011 @@ -23,7 +23,7 @@ import org.apache.commons.math.exception import org.apache.commons.math.exception.MathIllegalStateException; import org.apache.commons.math.linear.Array2DRowRealMatrix; import org.apache.commons.math.linear.RealMatrixPreservingVisitor; -import org.apache.commons.math.ode.FirstOrderDifferentialEquations; +import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations; import org.apache.commons.math.ode.sampling.NordsieckStepInterpolator; import org.apache.commons.math.ode.sampling.StepHandler; import org.apache.commons.math.util.FastMath; @@ -206,24 +206,29 @@ public class AdamsMoultonIntegrator exte /** {@inheritDoc} */ @Override - public double integrate(final FirstOrderDifferentialEquations equations, - final double t0, final double[] y0, - final double t, final double[] y) + public double integrate(final ExpandableFirstOrderDifferentialEquations equations, + final double t0, final double[] z0, + final double t, final double[] z) throws MathIllegalStateException, MathIllegalArgumentException { - final int n = y0.length; - sanityChecks(equations, t0, y0, t, y); + sanityChecks(equations, t0, z0, t, z); setEquations(equations); resetEvaluations(); final boolean forward = t > t0; // initialize working arrays + final int totalDim = equations.getDimension(); + final int mainDim = equations.getMainSetDimension(); + final double[] y0 = new double[totalDim]; + final double[] y = new double[totalDim]; + System.arraycopy(z0, 0, y0, 0, mainDim); + System.arraycopy(equations.getCurrentAdditionalStates(), 0, y0, mainDim, totalDim - mainDim); if (y != y0) { - System.arraycopy(y0, 0, y, 0, n); + System.arraycopy(y0, 0, y, 0, totalDim); } - final double[] yDot = new double[y0.length]; - final double[] yTmp = new double[y0.length]; - final double[] predictedScaled = new double[y0.length]; + final double[] yDot = new double[totalDim]; + final double[] yTmp = new double[totalDim]; + final double[] predictedScaled = new double[totalDim]; Array2DRowRealMatrix nordsieckTmp = null; // set up two interpolators sharing the integrator arrays @@ -290,7 +295,7 @@ public class AdamsMoultonIntegrator exte updateHighOrderDerivativesPhase2(predictedScaled, correctedScaled, nordsieckTmp); // discrete events handling - System.arraycopy(yTmp, 0, y, 0, n); + System.arraycopy(yTmp, 0, y, 0, totalDim); interpolator.reinitialize(stepEnd, stepSize, correctedScaled, nordsieckTmp); interpolator.storeTime(stepStart); interpolator.shift(); @@ -330,6 +335,10 @@ public class AdamsMoultonIntegrator exte } while (!isLastStep); + // dispatch result between main and additional states + System.arraycopy(y, 0, z, 0, z.length); + equations.setCurrentAdditionalState(y); + final double stopTime = stepStart; stepStart = Double.NaN; stepSize = Double.NaN; Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java?rev=1175409&r1=1175408&r2=1175409&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java Sun Sep 25 15:04:39 2011 @@ -23,7 +23,7 @@ import org.apache.commons.math.exception import org.apache.commons.math.exception.NumberIsTooSmallException; import org.apache.commons.math.exception.util.LocalizedFormats; import org.apache.commons.math.ode.AbstractIntegrator; -import org.apache.commons.math.ode.ExtendedFirstOrderDifferentialEquations; +import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations; import org.apache.commons.math.ode.FirstOrderDifferentialEquations; import org.apache.commons.math.util.FastMath; @@ -42,12 +42,11 @@ import org.apache.commons.math.util.Fast * component. The user can also use only two scalar values absTol and * relTol which will be used for all components. * </p> - * - * <p>If the Ordinary Differential Equations is an {@link ExtendedFirstOrderDifferentialEquations - * extended ODE} rather than a {@link FirstOrderDifferentialEquations basic ODE}, - * then <em>only</em> the {@link ExtendedFirstOrderDifferentialEquations#getMainSetDimension() - * main set} part of the state vector is used for stepsize control, not the complete - * state vector. + * <p> + * If the Ordinary Differential Equations is an {@link ExpandableFirstOrderDifferentialEquations + * extended ODE} rather than a {@link FirstOrderDifferentialEquations basic ODE}, then + * <em>only</em> the {@link ExpandableFirstOrderDifferentialEquations#getMainSet() main part} + * of the state vector is used for stepsize control, not the complete state vector. * </p> * * <p>If the estimated error for ym+1 is such that @@ -224,18 +223,14 @@ public abstract class AdaptiveStepsizeIn * @exception NumberIsTooSmallException if integration span is too small */ @Override - protected void sanityChecks(final FirstOrderDifferentialEquations equations, + protected void sanityChecks(final ExpandableFirstOrderDifferentialEquations equations, final double t0, final double[] y0, final double t, final double[] y) throws DimensionMismatchException, NumberIsTooSmallException { super.sanityChecks(equations, t0, y0, t, y); - if (equations instanceof ExtendedFirstOrderDifferentialEquations) { - mainSetDimension = ((ExtendedFirstOrderDifferentialEquations) equations).getMainSetDimension(); - } else { - mainSetDimension = equations.getDimension(); - } + mainSetDimension = equations.getMainSetDimension(); if ((vecAbsoluteTolerance != null) && (vecAbsoluteTolerance.length != mainSetDimension)) { throw new DimensionMismatchException(mainSetDimension, vecAbsoluteTolerance.length); @@ -248,7 +243,6 @@ public abstract class AdaptiveStepsizeIn } /** Initialize the integration step. - * @param equations differential equations set * @param forward forward integration indicator * @param order order of the method * @param scale scaling vector for the state vector (can be shorter than state vector) @@ -259,8 +253,7 @@ public abstract class AdaptiveStepsizeIn * @param yDot1 work array for the first time derivative of y1 * @return first integration step */ - public double initializeStep(final FirstOrderDifferentialEquations equations, - final boolean forward, final int order, final double[] scale, + public double initializeStep(final boolean forward, final int order, final double[] scale, final double t0, final double[] y0, final double[] yDot0, final double[] y1, final double[] yDot1) { @@ -356,7 +349,7 @@ public abstract class AdaptiveStepsizeIn } /** {@inheritDoc} */ - public abstract double integrate (FirstOrderDifferentialEquations equations, + public abstract double integrate (ExpandableFirstOrderDifferentialEquations equations, double t0, double[] y0, double t, double[] y) throws MathIllegalStateException, MathIllegalArgumentException; Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java?rev=1175409&r1=1175408&r2=1175409&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java Sun Sep 25 15:04:39 2011 @@ -19,7 +19,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.ode.FirstOrderDifferentialEquations; +import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations; import org.apache.commons.math.ode.sampling.StepHandler; import org.apache.commons.math.util.FastMath; @@ -189,24 +189,30 @@ public abstract class EmbeddedRungeKutta /** {@inheritDoc} */ @Override - public double integrate(final FirstOrderDifferentialEquations equations, - final double t0, final double[] y0, - final double t, final double[] y) + public double integrate(final ExpandableFirstOrderDifferentialEquations equations, + final double t0, final double[] z0, + final double t, final double[] z) throws MathIllegalStateException, MathIllegalArgumentException { - sanityChecks(equations, t0, y0, t, y); + sanityChecks(equations, t0, z0, t, z); setEquations(equations); resetEvaluations(); final boolean forward = t > t0; // create some internal working arrays + final int totalDim = equations.getDimension(); + final int mainDim = equations.getMainSetDimension(); + final double[] y0 = new double[totalDim]; + final double[] y = new double[totalDim]; + System.arraycopy(z0, 0, y0, 0, mainDim); + System.arraycopy(equations.getCurrentAdditionalStates(), 0, y0, mainDim, totalDim - mainDim); final int stages = c.length + 1; if (y != y0) { - System.arraycopy(y0, 0, y, 0, y0.length); + System.arraycopy(y0, 0, y, 0, totalDim); } - final double[][] yDotK = new double[stages][y0.length]; - final double[] yTmp = new double[y0.length]; - final double[] yDotTmp = new double[y0.length]; + final double[][] yDotK = new double[stages][totalDim]; + final double[] yTmp = new double[totalDim]; + final double[] yDotTmp = new double[totalDim]; // set up an interpolator sharing the integrator arrays final RungeKuttaStepInterpolator interpolator = (RungeKuttaStepInterpolator) prototype.copy(); @@ -248,7 +254,7 @@ public abstract class EmbeddedRungeKutta scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * FastMath.abs(y[i]); } } - hNew = initializeStep(equations, forward, getOrder(), scale, + hNew = initializeStep(forward, getOrder(), scale, stepStart, y, yDotK[0], yTmp, yDotK[1]); firstTime = false; } @@ -325,6 +331,10 @@ public abstract class EmbeddedRungeKutta } 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; Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java?rev=1175409&r1=1175408&r2=1175409&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java Sun Sep 25 15:04:39 2011 @@ -20,7 +20,7 @@ package org.apache.commons.math.ode.nons import org.apache.commons.math.analysis.solvers.UnivariateRealSolver; import org.apache.commons.math.exception.MathIllegalArgumentException; import org.apache.commons.math.exception.MathIllegalStateException; -import org.apache.commons.math.ode.FirstOrderDifferentialEquations; +import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations; import org.apache.commons.math.ode.events.EventHandler; import org.apache.commons.math.ode.sampling.AbstractStepInterpolator; import org.apache.commons.math.ode.sampling.StepHandler; @@ -541,26 +541,32 @@ public class GraggBulirschStoerIntegrato /** {@inheritDoc} */ @Override - public double integrate(final FirstOrderDifferentialEquations equations, - final double t0, final double[] y0, final double t, final double[] y) + public double integrate(final ExpandableFirstOrderDifferentialEquations equations, + final double t0, final double[] z0, final double t, final double[] z) throws MathIllegalStateException, MathIllegalArgumentException { - sanityChecks(equations, t0, y0, t, y); + sanityChecks(equations, t0, z0, t, z); setEquations(equations); resetEvaluations(); final boolean forward = t > t0; // create some internal working arrays - final double[] yDot0 = new double[y0.length]; - final double[] y1 = new double[y0.length]; - final double[] yTmp = new double[y0.length]; - final double[] yTmpDot = new double[y0.length]; + final int totalDim = equations.getDimension(); + final int mainDim = equations.getMainSetDimension(); + final double[] y0 = new double[totalDim]; + final double[] y = new double[totalDim]; + System.arraycopy(z0, 0, y0, 0, mainDim); + System.arraycopy(equations.getCurrentAdditionalStates(), 0, y0, mainDim, totalDim - mainDim); + final double[] yDot0 = new double[totalDim]; + final double[] y1 = new double[totalDim]; + final double[] yTmp = new double[totalDim]; + final double[] yTmpDot = new double[totalDim]; final double[][] diagonal = new double[sequence.length-1][]; final double[][] y1Diag = new double[sequence.length-1][]; for (int k = 0; k < sequence.length-1; ++k) { - diagonal[k] = new double[y0.length]; - y1Diag[k] = new double[y0.length]; + diagonal[k] = new double[totalDim]; + y1Diag[k] = new double[totalDim]; } final double[][][] fk = new double[sequence.length][][]; @@ -631,8 +637,7 @@ public class GraggBulirschStoerIntegrato } if (firstTime) { - hNew = initializeStep(equations, forward, - 2 * targetIter + 1, scale, + hNew = initializeStep(forward, 2 * targetIter + 1, scale, stepStart, y, yDot0, yTmp, yTmpDot); } Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java?rev=1175409&r1=1175408&r2=1175409&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java Sun Sep 25 15:04:39 2011 @@ -21,7 +21,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.ode.AbstractIntegrator; -import org.apache.commons.math.ode.FirstOrderDifferentialEquations; +import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations; import org.apache.commons.math.ode.sampling.StepHandler; import org.apache.commons.math.util.FastMath; @@ -90,17 +90,23 @@ public abstract class RungeKuttaIntegrat } /** {@inheritDoc} */ - public double integrate(final FirstOrderDifferentialEquations equations, - final double t0, final double[] y0, - final double t, final double[] y) + public double integrate(final ExpandableFirstOrderDifferentialEquations equations, + final double t0, final double[] z0, + final double t, final double[] z) throws MathIllegalStateException, MathIllegalArgumentException { - sanityChecks(equations, t0, y0, t, y); + sanityChecks(equations, t0, z0, t, z); setEquations(equations); resetEvaluations(); final boolean forward = t > t0; // create some internal working arrays + final int totalDim = equations.getDimension(); + final int mainDim = equations.getMainSetDimension(); + final double[] y0 = new double[totalDim]; + final double[] y = new double[totalDim]; + System.arraycopy(z0, 0, y0, 0, mainDim); + System.arraycopy(equations.getCurrentAdditionalStates(), 0, y0, mainDim, totalDim - mainDim); final int stages = c.length + 1; if (y != y0) { System.arraycopy(y0, 0, y, 0, y0.length); @@ -179,6 +185,10 @@ public abstract class RungeKuttaIntegrat } while (!isLastStep); + // dispatch result between main and additional states + System.arraycopy(y, 0, z, 0, z.length); + equations.setCurrentAdditionalState(y); + final double stopTime = stepStart; stepStart = Double.NaN; stepSize = Double.NaN; Modified: commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties?rev=1175409&r1=1175408&r2=1175409&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties (original) +++ commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties Sun Sep 25 15:04:39 2011 @@ -279,6 +279,9 @@ UNABLE_TO_PERFORM_QR_DECOMPOSITION_ON_JA UNABLE_TO_SOLVE_SINGULAR_PROBLEM = r\u00e9solution impossible : probl\u00e8me singulier UNBOUNDED_SOLUTION = solution non born\u00e9e UNKNOWN_MODE = mode {0} inconnu, modes connus : {1} ({2}), {3} ({4}), {5} ({6}), {7} ({8}), {9} ({10}) et {11} ({12}) +UNKNOWN_ADDITIONAL_EQUATION = \u00e9quation additionnelle inconnue +UNKNOWN_PARAMETER = param\u00e8tre {0} inconnu +UNMATCHED_ODE_IN_EXTENDED_SET = l''\u00e9quation diff\u00e9rentielle ne correspond pas \u00e0 l''\u00e9quation principale du jeu \u00e9tendu CANNOT_PARSE_AS_TYPE = cha\u00eene {0} non analysable (\u00e0 partir de la position {1}) en un objet de type {2} CANNOT_PARSE = cha\u00eene {0} non analysable (\u00e0 partir de la position {1}) UNPARSEABLE_3D_VECTOR = vecteur 3D non analysable : "{0}" 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=1175409&r1=1175408&r2=1175409&view=diff ============================================================================== --- commons/proper/math/trunk/src/site/xdoc/changes.xml (original) +++ commons/proper/math/trunk/src/site/xdoc/changes.xml Sun Sep 25 15:04:39 2011 @@ -52,6 +52,12 @@ The <action> type attribute can be add,u If the output is not quite correct, check for invisible trailing spaces! --> <release version="3.0" date="TBD" description="TBD"> + <action dev="luc" type="fix" issue="MATH-380" due-to="Pascal Parraud"> + Completely revamped the computation of Jacobians in ODE. This computation is now + included in the mainstream class hierarchy, not in a separate package anymore, + and it allows adding other types of equations to a main ODE, not only variational + equations for Jacobians computation. + </action> <action dev="gregs" type="update" issue="MATH-607"> SimpleRegression implements UpdatingMultipleLinearRegression interface. </action> Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/JacobianMatricesTest.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/JacobianMatricesTest.java?rev=1175409&view=auto ============================================================================== --- commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/JacobianMatricesTest.java (added) +++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/JacobianMatricesTest.java Sun Sep 25 15:04:39 2011 @@ -0,0 +1,598 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math.ode; + +import org.apache.commons.math.ode.nonstiff.DormandPrince54Integrator; +import org.apache.commons.math.stat.descriptive.SummaryStatistics; +import org.apache.commons.math.util.FastMath; +import org.junit.Assert; +import org.junit.Test; + +public class JacobianMatricesTest { + + @Test + public void testLowAccuracyExternalDifferentiation() { + // this test does not really test FirstOrderIntegratorWithJacobians, + // it only shows that WITHOUT this class, attempting to recover + // the jacobians from external differentiation on simple integration + // results with low accuracy gives very poor results. In fact, + // the curves dy/dp = g(b) when b varies from 2.88 to 3.08 are + // essentially noise. + // This test is taken from Hairer, Norsett and Wanner book + // Solving Ordinary Differential Equations I (Nonstiff problems), + // the curves dy/dp = g(b) are in figure 6.5 + FirstOrderIntegrator integ = + new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 }); + double hP = 1.0e-12; + SummaryStatistics residualsP0 = new SummaryStatistics(); + SummaryStatistics residualsP1 = new SummaryStatistics(); + for (double b = 2.88; b < 3.08; b += 0.001) { + Brusselator brusselator = new Brusselator(b); + double[] y = { 1.3, b }; + integ.integrate(brusselator, 0, y, 20.0, y); + double[] yP = { 1.3, b + hP }; + integ.integrate(brusselator, 0, yP, 20.0, yP); + residualsP0.addValue((yP[0] - y[0]) / hP - brusselator.dYdP0()); + residualsP1.addValue((yP[1] - y[1]) / hP - brusselator.dYdP1()); + } + Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) > 500); + Assert.assertTrue(residualsP0.getStandardDeviation() > 30); + Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) > 700); + Assert.assertTrue(residualsP1.getStandardDeviation() > 40); + } + + @Test + public void testHighAccuracyExternalDifferentiation() { + FirstOrderIntegrator integ = + new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 1.0e-10, 1.0e-10 }); + double hP = 1.0e-12; + SummaryStatistics residualsP0 = new SummaryStatistics(); + SummaryStatistics residualsP1 = new SummaryStatistics(); + for (double b = 2.88; b < 3.08; b += 0.001) { + ParamBrusselator brusselator = new ParamBrusselator(b); + double[] y = { 1.3, b }; + integ.integrate(brusselator, 0, y, 20.0, y); + double[] yP = { 1.3, b + hP }; + brusselator.setParameter("b", b + hP); + integ.integrate(brusselator, 0, yP, 20.0, yP); + residualsP0.addValue((yP[0] - y[0]) / hP - brusselator.dYdP0()); + residualsP1.addValue((yP[1] - y[1]) / hP - brusselator.dYdP1()); + } + Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) > 0.02); + Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.03); + Assert.assertTrue(residualsP0.getStandardDeviation() > 0.003); + Assert.assertTrue(residualsP0.getStandardDeviation() < 0.004); + Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) > 0.04); + Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05); + Assert.assertTrue(residualsP1.getStandardDeviation() > 0.007); + Assert.assertTrue(residualsP1.getStandardDeviation() < 0.008); + } + + @Test + public void testInternalDifferentiation() { + ExpandableFirstOrderIntegrator integ = + new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 }); + double hP = 1.0e-12; + double hY = 1.0e-12; + SummaryStatistics residualsP0 = new SummaryStatistics(); + SummaryStatistics residualsP1 = new SummaryStatistics(); + for (double b = 2.88; b < 3.08; b += 0.001) { + ParamBrusselator brusselator = new ParamBrusselator(b); + brusselator.setParameter(ParamBrusselator.B, b); + double[] z = { 1.3, b }; + double[][] dZdZ0 = new double[2][2]; + double[] dZdP = new double[2]; + ExpandableFirstOrderDifferentialEquations efode = new ExpandableFirstOrderDifferentialEquations(brusselator); + JacobianMatrices jacob = new JacobianMatrices(efode, brusselator); + jacob.setParameterizedODE(brusselator); + jacob.selectParameters(ParamBrusselator.B); + jacob.setMainStateSteps(new double[] { hY, hY }); + jacob.setParameterStep(ParamBrusselator.B, hP); + jacob.setInitialParameterJacobian(ParamBrusselator.B, new double[] { 0.0, 1.0 }); + integ.setMaxEvaluations(5000); + integ.integrate(efode, 0, z, 20.0, z); + jacob.getCurrentMainSetJacobian(dZdZ0); + jacob.getCurrentParameterJacobian(ParamBrusselator.B, dZdP); +// Assert.assertEquals(5000, integ.getMaxEvaluations()); +// Assert.assertTrue(integ.getEvaluations() > 1500); +// Assert.assertTrue(integ.getEvaluations() < 2100); +// Assert.assertEquals(4 * integ.getEvaluations(), integ.getEvaluations()); + residualsP0.addValue(dZdP[0] - brusselator.dYdP0()); + residualsP1.addValue(dZdP[1] - brusselator.dYdP1()); + } + Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.02); + Assert.assertTrue(residualsP0.getStandardDeviation() < 0.003); + Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05); + Assert.assertTrue(residualsP1.getStandardDeviation() < 0.01); + } + + @Test + public void testAnalyticalDifferentiation() { + ExpandableFirstOrderIntegrator integ = + new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 }); + SummaryStatistics residualsP0 = new SummaryStatistics(); + SummaryStatistics residualsP1 = new SummaryStatistics(); + for (double b = 2.88; b < 3.08; b += 0.001) { + Brusselator brusselator = new Brusselator(b); + double[] z = { 1.3, b }; + double[][] dZdZ0 = new double[2][2]; + double[] dZdP = new double[2]; + ExpandableFirstOrderDifferentialEquations efode = new ExpandableFirstOrderDifferentialEquations(brusselator); + JacobianMatrices jacob = new JacobianMatrices(efode, brusselator); + jacob.setParameterJacobianProvider(brusselator); + jacob.selectParameters(Brusselator.B); + jacob.setInitialParameterJacobian(Brusselator.B, new double[] { 0.0, 1.0 }); + integ.setMaxEvaluations(5000); + integ.integrate(efode, 0, z, 20.0, z); + jacob.getCurrentMainSetJacobian(dZdZ0); + jacob.getCurrentParameterJacobian(Brusselator.B, dZdP); +// Assert.assertEquals(5000, integ.getMaxEvaluations()); +// Assert.assertTrue(integ.getEvaluations() > 350); +// Assert.assertTrue(integ.getEvaluations() < 510); + residualsP0.addValue(dZdP[0] - brusselator.dYdP0()); + residualsP1.addValue(dZdP[1] - brusselator.dYdP1()); + } + Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.014); + Assert.assertTrue(residualsP0.getStandardDeviation() < 0.003); + Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05); + Assert.assertTrue(residualsP1.getStandardDeviation() < 0.01); + } + + @Test + public void testFinalResult() { + + ExpandableFirstOrderIntegrator integ = + new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 1.0e-10, 1.0e-10 }); + double[] y = new double[] { 0.0, 1.0 }; + Circle circle = new Circle(y, 1.0, 1.0, 0.1); + + ExpandableFirstOrderDifferentialEquations efode = new ExpandableFirstOrderDifferentialEquations(circle); + JacobianMatrices jacob = new JacobianMatrices(efode, circle); + jacob.setParameterJacobianProvider(circle); + jacob.selectParameters(Circle.CX, Circle.CY, Circle.OMEGA); + jacob.setInitialMainStateJacobian(circle.exactDyDy0(0)); + jacob.setInitialParameterJacobian(Circle.CX, circle.exactDyDcx(0)); + jacob.setInitialParameterJacobian(Circle.CY, circle.exactDyDcy(0)); + jacob.setInitialParameterJacobian(Circle.OMEGA, circle.exactDyDom(0)); + integ.setMaxEvaluations(5000); + + double t = 18 * FastMath.PI; + integ.integrate(efode, 0, y, t, y); + for (int i = 0; i < y.length; ++i) { + Assert.assertEquals(circle.exactY(t)[i], y[i], 1.0e-9); + } + + double[][] dydy0 = new double[2][2]; + jacob.getCurrentMainSetJacobian(dydy0); + for (int i = 0; i < dydy0.length; ++i) { + for (int j = 0; j < dydy0[i].length; ++j) { + Assert.assertEquals(circle.exactDyDy0(t)[i][j], dydy0[i][j], 1.0e-9); + } + } + double[] dydcx = new double[2]; + jacob.getCurrentParameterJacobian(Circle.CX, dydcx); + for (int i = 0; i < dydcx.length; ++i) { + Assert.assertEquals(circle.exactDyDcx(t)[i], dydcx[i], 1.0e-7); + } + double[] dydcy = new double[2]; + jacob.getCurrentParameterJacobian(Circle.CY, dydcy); + for (int i = 0; i < dydcy.length; ++i) { + Assert.assertEquals(circle.exactDyDcy(t)[i], dydcy[i], 1.0e-7); + } + double[] dydom = new double[2]; + jacob.getCurrentParameterJacobian(Circle.OMEGA, dydom); + for (int i = 0; i < dydom.length; ++i) { + Assert.assertEquals(circle.exactDyDom(t)[i], dydom[i], 1.0e-7); + } + } + + @Test + public void testParameterizable() { + + ExpandableFirstOrderIntegrator integ = + new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 1.0e-10, 1.0e-10 }); + double[] y = new double[] { 0.0, 1.0 }; + ParameterizedCircle pcircle = new ParameterizedCircle(y, 1.0, 1.0, 0.1); +// pcircle.setParameter(ParameterizedCircle.CX, 1.0); +// pcircle.setParameter(ParameterizedCircle.CY, 1.0); +// pcircle.setParameter(ParameterizedCircle.OMEGA, 0.1); + + ExpandableFirstOrderDifferentialEquations efode = new ExpandableFirstOrderDifferentialEquations(pcircle); + + double hP = 1.0e-12; + double hY = 1.0e-12; + + JacobianMatrices jacob = new JacobianMatrices(efode, pcircle); + jacob.setParameterJacobianProvider(pcircle); + jacob.setParameterizedODE(pcircle); + jacob.selectParameters(Circle.CX, Circle.OMEGA); + jacob.setMainStateSteps(new double[] { hY, hY }); + jacob.setParameterStep(Circle.OMEGA, hP); + jacob.setInitialMainStateJacobian(pcircle.exactDyDy0(0)); + jacob.setInitialParameterJacobian(Circle.CX, pcircle.exactDyDcx(0)); +// jacob.setInitialParameterJacobian(Circle.CY, circle.exactDyDcy(0)); + jacob.setInitialParameterJacobian(Circle.OMEGA, pcircle.exactDyDom(0)); + integ.setMaxEvaluations(50000); + + double t = 18 * FastMath.PI; + integ.integrate(efode, 0, y, t, y); + for (int i = 0; i < y.length; ++i) { + Assert.assertEquals(pcircle.exactY(t)[i], y[i], 1.0e-9); + } + + double[][] dydy0 = new double[2][2]; + jacob.getCurrentMainSetJacobian(dydy0); + for (int i = 0; i < dydy0.length; ++i) { + for (int j = 0; j < dydy0[i].length; ++j) { + Assert.assertEquals(pcircle.exactDyDy0(t)[i][j], dydy0[i][j], 5.0e-4); + } + } + + double[] dydp0 = new double[2]; + jacob.getCurrentParameterJacobian(Circle.CX, dydp0); + for (int i = 0; i < dydp0.length; ++i) { + Assert.assertEquals(pcircle.exactDyDcx(t)[i], dydp0[i], 5.0e-4); + } + + double[] dydp1 = new double[2]; + jacob.getCurrentParameterJacobian(Circle.OMEGA, dydp1); + for (int i = 0; i < dydp1.length; ++i) { + Assert.assertEquals(pcircle.exactDyDom(t)[i], dydp1[i], 1.0e-2); + } + } + + private static class Brusselator extends AbstractParameterizable + implements MainStateJacobianProvider, ParameterJacobianProvider { + + public static final String B = "b"; + + private double b; + + public Brusselator(double b) { + super(B); + this.b = b; + } + + public int getDimension() { + return 2; + } + + public void computeDerivatives(double t, double[] y, double[] yDot) { + double prod = y[0] * y[0] * y[1]; + yDot[0] = 1 + prod - (b + 1) * y[0]; + yDot[1] = b * y[0] - prod; + } + + public void computeMainStateJacobian(double t, double[] y, double[] yDot, + double[][] dFdY) { + double p = 2 * y[0] * y[1]; + double y02 = y[0] * y[0]; + dFdY[0][0] = p - (1 + b); + dFdY[0][1] = y02; + dFdY[1][0] = b - p; + dFdY[1][1] = -y02; + } + + public void computeParameterJacobian(double t, double[] y, double[] yDot, + String paramName, double[] dFdP) { + complainIfNotSupported(paramName); + dFdP[0] = -y[0]; + dFdP[1] = y[0]; + } + + public double dYdP0() { + return -1088.232716447743 + (1050.775747149553 + (-339.012934631828 + 36.52917025056327 * b) * b) * b; + } + + public double dYdP1() { + return 1502.824469929139 + (-1438.6974831849952 + (460.959476642384 - 49.43847385647082 * b) * b) * b; + } + + } + + private static class ParamBrusselator extends AbstractParameterizable + implements FirstOrderDifferentialEquations, ParameterizedODE { + + public static final String B = "b"; + + private double b; + + public ParamBrusselator(double b) { + super(B); + this.b = b; + } + + public int getDimension() { + return 2; + } + + /** {@inheritDoc} */ + public double getParameter(final String name) + throws IllegalArgumentException { + complainIfNotSupported(name); + return b; + } + + /** {@inheritDoc} */ + public void setParameter(final String name, final double value) + throws IllegalArgumentException { + complainIfNotSupported(name); + b = value; + } + + public void computeDerivatives(double t, double[] y, double[] yDot) { + double prod = y[0] * y[0] * y[1]; + yDot[0] = 1 + prod - (b + 1) * y[0]; + yDot[1] = b * y[0] - prod; + } + + public double dYdP0() { + return -1088.232716447743 + (1050.775747149553 + (-339.012934631828 + 36.52917025056327 * b) * b) * b; + } + + public double dYdP1() { + return 1502.824469929139 + (-1438.6974831849952 + (460.959476642384 - 49.43847385647082 * b) * b) * b; + } + + } + + /** ODE representing a point moving on a circle with provided center and angular rate. */ + private static class Circle extends AbstractParameterizable + implements MainStateJacobianProvider, ParameterJacobianProvider { + + public static final String CX = "cx"; + public static final String CY = "cy"; + public static final String OMEGA = "omega"; + + private final double[] y0; + private double cx; + private double cy; + private double omega; + + public Circle(double[] y0, double cx, double cy, double omega) { + super(CX,CY,OMEGA); + this.y0 = y0.clone(); + this.cx = cx; + this.cy = cy; + this.omega = omega; + } + + public int getDimension() { + return 2; + } + + public void computeDerivatives(double t, double[] y, double[] yDot) { + yDot[0] = omega * (cy - y[1]); + yDot[1] = omega * (y[0] - cx); + } + + public void computeMainStateJacobian(double t, double[] y, + double[] yDot, double[][] dFdY) { + dFdY[0][0] = 0; + dFdY[0][1] = -omega; + dFdY[1][0] = omega; + dFdY[1][1] = 0; + } + + public void computeParameterJacobian(double t, double[] y, double[] yDot, + String paramName, double[] dFdP) { + complainIfNotSupported(paramName); + if (paramName.equals(CX)) { + dFdP[0] = 0; + dFdP[1] = -omega; + } else if (paramName.equals(CY)) { + dFdP[0] = omega; + dFdP[1] = 0; + } else { + dFdP[0] = cy - y[1]; + dFdP[1] = y[0] - cx; + } + } + + public double[] exactY(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + double dx0 = y0[0] - cx; + double dy0 = y0[1] - cy; + return new double[] { + cx + cos * dx0 - sin * dy0, + cy + sin * dx0 + cos * dy0 + }; + } + + public double[][] exactDyDy0(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + return new double[][] { + { cos, -sin }, + { sin, cos } + }; + } + + public double[] exactDyDcx(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + return new double[] {1 - cos, -sin}; + } + + public double[] exactDyDcy(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + return new double[] {sin, 1 - cos}; + } + + public double[] exactDyDom(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + double dx0 = y0[0] - cx; + double dy0 = y0[1] - cy; + return new double[] { -t * (sin * dx0 + cos * dy0) , t * (cos * dx0 - sin * dy0) }; + } + + public double[][] exactDyDp(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + double dx0 = y0[0] - cx; + double dy0 = y0[1] - cy; + return new double[][] { + { 1 - cos, sin, -t * (sin * dx0 + cos * dy0) }, + { -sin, 1 - cos, t * (cos * dx0 - sin * dy0) } + }; + } + + public double[] exactYDot(double t) { + double oCos = omega * FastMath.cos(omega * t); + double oSin = omega * FastMath.sin(omega * t); + double dx0 = y0[0] - cx; + double dy0 = y0[1] - cy; + return new double[] { + -oSin * dx0 - oCos * dy0, + oCos * dx0 - oSin * dy0 + }; + } + + public double[][] exactDyDy0Dot(double t) { + double oCos = omega * FastMath.cos(omega * t); + double oSin = omega * FastMath.sin(omega * t); + return new double[][] { + { -oSin, -oCos }, + { oCos, -oSin } + }; + } + + public double[][] exactDyDpDot(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + double oCos = omega * cos; + double oSin = omega * sin; + double dx0 = y0[0] - cx; + double dy0 = y0[1] - cy; + return new double[][] { + { oSin, oCos, -sin * dx0 - cos * dy0 - t * ( oCos * dx0 - oSin * dy0) }, + { -oCos, oSin, cos * dx0 - sin * dy0 + t * (-oSin * dx0 - oCos * dy0) } + }; + } + + } + + /** ODE representing a point moving on a circle with provided center and angular rate. */ + private static class ParameterizedCircle extends AbstractParameterizable + implements FirstOrderDifferentialEquations, ParameterJacobianProvider, ParameterizedODE { + + public static final String CX = "cx"; + public static final String CY = "cy"; + public static final String OMEGA = "omega"; + + private final double[] y0; + private double cx; + private double cy; + private double omega; + + public ParameterizedCircle(double[] y0, double cx, double cy, double omega) { + super(CX,CY,OMEGA); + this.y0 = y0.clone(); + this.cx = cx; + this.cy = cy; + this.omega = omega; + } + + public int getDimension() { + return 2; + } + + public void computeDerivatives(double t, double[] y, double[] yDot) { + yDot[0] = omega * (cy - y[1]); + yDot[1] = omega * (y[0] - cx); + } + + public void computeParameterJacobian(double t, double[] y, double[] yDot, + String paramName, double[] dFdP) + throws IllegalArgumentException { + if (paramName.equals(CX)) { + dFdP[0] = 0; + dFdP[1] = -omega; + } else { + throw new IllegalArgumentException(); + } + } + + public double getParameter(final String name) + throws IllegalArgumentException { + if (name.equals(CY)) { + return cy; + } else if (name.equals(OMEGA)) { + return omega; + } else { + throw new IllegalArgumentException(); + } + } + + public void setParameter(final String name, final double value) + throws IllegalArgumentException { + if (name.equals(CY)) { + cy = value; + } else if (name.equals(OMEGA)) { + omega = value; + } else { + throw new IllegalArgumentException(); + } + } + + public double[] exactY(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + double dx0 = y0[0] - cx; + double dy0 = y0[1] - cy; + return new double[] { + cx + cos * dx0 - sin * dy0, + cy + sin * dx0 + cos * dy0 + }; + } + + public double[][] exactDyDy0(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + return new double[][] { + { cos, -sin }, + { sin, cos } + }; + } + + public double[] exactDyDcx(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + return new double[] {1 - cos, -sin}; + } + + public double[] exactDyDcy(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + return new double[] {sin, 1 - cos}; + } + + public double[] exactDyDom(double t) { + double cos = FastMath.cos(omega * t); + double sin = FastMath.sin(omega * t); + double dx0 = y0[0] - cx; + double dy0 = y0[1] - cy; + return new double[] { -t * (sin * dx0 + cos * dy0) , t * (cos * dx0 - sin * dy0) }; + } + + } + +} Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/JacobianMatricesTest.java ------------------------------------------------------------------------------ svn:eol-style = native Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/JacobianMatricesTest.java ------------------------------------------------------------------------------ svn:keywords = "Author Date Id Revision" Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java?rev=1175409&r1=1175408&r2=1175409&view=diff ============================================================================== --- commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java (original) +++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java Sun Sep 25 15:04:39 2011 @@ -256,8 +256,7 @@ public class GraggBulirschStoerIntegrato } @Test - public void testUnstableDerivative() - { + public void testUnstableDerivative() { final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0); FirstOrderIntegrator integ = new GraggBulirschStoerIntegrator(0.1, 10, 1.0e-12, 0.0);
