Added a feature allowing error estimation to be computed only on a subset of
Ordinary Differential Equations, considered as the main set, the remaining equations being considered only as an extension set that should not influence the ODE integration algorithm JIRA: MATH-388 git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@967007 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
730e17083e
commit
92f005ad13
|
|
@ -0,0 +1,66 @@
|
|||
/*
|
||||
* 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 a first order differential equations set
|
||||
* with a main set of equations and an extension set.
|
||||
*
|
||||
* <p>
|
||||
* This interface is a simple extension on the {@link
|
||||
* FirstOrderDifferentialEquations} that allows to identify which part
|
||||
* of a complete set of differential equations correspond to the main
|
||||
* set and which part correspond to the extension set.
|
||||
* </p>
|
||||
* <p>
|
||||
* One typical use case is the computation of Jacobians. The main
|
||||
* set of equations correspond to the raw ode, and we add to this set
|
||||
* another bunch of equations which represent the jacobians of the
|
||||
* main set. In that case, we want the integrator to use <em>only</em>
|
||||
* the main set to estimate the errors and hence the step sizes. It should
|
||||
* <em>not</em> use the additional equations in this computation. If the
|
||||
* complete ode implements this interface, the {@link FirstOrderIntegrator
|
||||
* integrator} will be able to know where the main set ends and where the
|
||||
* extended set begins.
|
||||
* </p>
|
||||
* <p>
|
||||
* We consider that the main set always corresponds to the first equations
|
||||
* and the extended set to the last equations.
|
||||
* </p>
|
||||
*
|
||||
* @see FirstOrderDifferentialEquations
|
||||
*
|
||||
* @version $Revision$ $Date$
|
||||
* @since 2.2
|
||||
*/
|
||||
|
||||
public interface ExtendedFirstOrderDifferentialEquations extends FirstOrderDifferentialEquations {
|
||||
|
||||
/** Return the dimension of the main set of equations.
|
||||
* <p>
|
||||
* The main set of equations represent 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 represent an extension set of equations which are considered
|
||||
* secondary.
|
||||
* </p>
|
||||
* @return dimension of the main set of equations, must be lesser than or
|
||||
* equal to the {@link #getDimension() total dimension}
|
||||
*/
|
||||
int getMainSetDimension();
|
||||
|
||||
}
|
||||
|
|
@ -378,7 +378,7 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
|
|||
}
|
||||
|
||||
/** Wrapper for differential equations, ensuring start evaluations are counted. */
|
||||
private class CountingDifferentialEquations implements FirstOrderDifferentialEquations {
|
||||
private class CountingDifferentialEquations implements ExtendedFirstOrderDifferentialEquations {
|
||||
|
||||
/** Dimension of the problem. */
|
||||
private final int dimension;
|
||||
|
|
@ -400,6 +400,11 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
|
|||
public int getDimension() {
|
||||
return dimension;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public int getMainSetDimension() {
|
||||
return mainSetDimension;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -27,6 +27,7 @@ import java.util.Collection;
|
|||
import org.apache.commons.math.MathRuntimeException;
|
||||
import org.apache.commons.math.MaxEvaluationsExceededException;
|
||||
import org.apache.commons.math.ode.DerivativeException;
|
||||
import org.apache.commons.math.ode.ExtendedFirstOrderDifferentialEquations;
|
||||
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
|
||||
import org.apache.commons.math.ode.FirstOrderIntegrator;
|
||||
import org.apache.commons.math.ode.IntegratorException;
|
||||
|
|
@ -354,7 +355,7 @@ public class FirstOrderIntegratorWithJacobians {
|
|||
}
|
||||
|
||||
/** Wrapper class used to map state and jacobians into compound state. */
|
||||
private class MappingWrapper implements FirstOrderDifferentialEquations {
|
||||
private class MappingWrapper implements ExtendedFirstOrderDifferentialEquations {
|
||||
|
||||
/** Current state. */
|
||||
private final double[] y;
|
||||
|
|
@ -388,6 +389,11 @@ public class FirstOrderIntegratorWithJacobians {
|
|||
return n * (1 + n + k);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public int getMainSetDimension() {
|
||||
return ode.getDimension();
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public void computeDerivatives(final double t, final double[] z, final double[] zDot)
|
||||
throws DerivativeException {
|
||||
|
|
@ -442,8 +448,7 @@ public class FirstOrderIntegratorWithJacobians {
|
|||
}
|
||||
|
||||
/** Wrapper class to compute jacobians by finite differences for ODE which do not compute them themselves. */
|
||||
private class FiniteDifferencesWrapper
|
||||
implements ODEWithJacobians {
|
||||
private class FiniteDifferencesWrapper implements ODEWithJacobians {
|
||||
|
||||
/** Raw ODE without jacobians computation. */
|
||||
private final ParameterizedODE ode;
|
||||
|
|
|
|||
|
|
@ -29,8 +29,7 @@ import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
|
|||
* @since 2.1
|
||||
*/
|
||||
|
||||
public interface ParameterizedODE
|
||||
extends FirstOrderDifferentialEquations {
|
||||
public interface ParameterizedODE extends FirstOrderDifferentialEquations {
|
||||
|
||||
/** Get the number of parameters.
|
||||
* @return number of parameters
|
||||
|
|
|
|||
|
|
@ -235,7 +235,7 @@ public class AdamsBashforthIntegrator extends AdamsIntegrator {
|
|||
|
||||
// evaluate error using the last term of the Taylor expansion
|
||||
error = 0;
|
||||
for (int i = 0; i < y0.length; ++i) {
|
||||
for (int i = 0; i < mainSetDimension; ++i) {
|
||||
final double yScale = Math.abs(y[i]);
|
||||
final double tol = (vecAbsoluteTolerance == null) ?
|
||||
(scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
|
||||
|
|
@ -243,7 +243,7 @@ public class AdamsBashforthIntegrator extends AdamsIntegrator {
|
|||
final double ratio = nordsieck.getEntry(lastRow, i) / tol;
|
||||
error += ratio * ratio;
|
||||
}
|
||||
error = Math.sqrt(error / y0.length);
|
||||
error = Math.sqrt(error / mainSetDimension);
|
||||
|
||||
if (error <= 1.0) {
|
||||
|
||||
|
|
|
|||
|
|
@ -419,7 +419,7 @@ public class AdamsMoultonIntegrator extends AdamsIntegrator {
|
|||
}
|
||||
|
||||
/**
|
||||
* End visiting te Nordsieck vector.
|
||||
* End visiting the Nordsieck vector.
|
||||
* <p>The correction is used to control stepsize. So its amplitude is
|
||||
* considered to be an error, which must be normalized according to
|
||||
* error control settings. If the normalized value is greater than 1,
|
||||
|
|
@ -432,15 +432,17 @@ public class AdamsMoultonIntegrator extends AdamsIntegrator {
|
|||
double error = 0;
|
||||
for (int i = 0; i < after.length; ++i) {
|
||||
after[i] += previous[i] + scaled[i];
|
||||
final double yScale = Math.max(Math.abs(previous[i]), Math.abs(after[i]));
|
||||
final double tol = (vecAbsoluteTolerance == null) ?
|
||||
(scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
|
||||
(vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yScale);
|
||||
final double ratio = (after[i] - before[i]) / tol;
|
||||
error += ratio * ratio;
|
||||
if (i < mainSetDimension) {
|
||||
final double yScale = Math.max(Math.abs(previous[i]), Math.abs(after[i]));
|
||||
final double tol = (vecAbsoluteTolerance == null) ?
|
||||
(scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
|
||||
(vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yScale);
|
||||
final double ratio = (after[i] - before[i]) / tol;
|
||||
error += ratio * ratio;
|
||||
}
|
||||
}
|
||||
|
||||
return Math.sqrt(error / after.length);
|
||||
return Math.sqrt(error / mainSetDimension);
|
||||
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -19,6 +19,7 @@ package org.apache.commons.math.ode.nonstiff;
|
|||
|
||||
import org.apache.commons.math.ode.AbstractIntegrator;
|
||||
import org.apache.commons.math.ode.DerivativeException;
|
||||
import org.apache.commons.math.ode.ExtendedFirstOrderDifferentialEquations;
|
||||
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
|
||||
import org.apache.commons.math.ode.IntegratorException;
|
||||
import org.apache.commons.math.util.LocalizedFormats;
|
||||
|
|
@ -36,14 +37,22 @@ import org.apache.commons.math.util.LocalizedFormats;
|
|||
* where absTol_i is the absolute tolerance for component i of the
|
||||
* state vector and relTol_i is the relative tolerance for the same
|
||||
* component. The user can also use only two scalar values absTol and
|
||||
* relTol which will be used for all components.</p>
|
||||
* 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>
|
||||
*
|
||||
* <p>If the estimated error for ym+1 is such that
|
||||
* <pre>
|
||||
* sqrt((sum (errEst_i / threshold_i)^2 ) / n) < 1
|
||||
* </pre>
|
||||
*
|
||||
* (where n is the state vector dimension) then the step is accepted,
|
||||
* (where n is the main set dimension) then the step is accepted,
|
||||
* otherwise the step is rejected and a new attempt is made with a new
|
||||
* stepsize.</p>
|
||||
*
|
||||
|
|
@ -76,6 +85,9 @@ public abstract class AdaptiveStepsizeIntegrator
|
|||
/** Maximal step. */
|
||||
private final double maxStep;
|
||||
|
||||
/** Main set dimension. */
|
||||
protected int mainSetDimension;
|
||||
|
||||
/** Build an integrator with the given stepsize bounds.
|
||||
* The default step handler does nothing.
|
||||
* @param name name of the method
|
||||
|
|
@ -171,14 +183,20 @@ public abstract class AdaptiveStepsizeIntegrator
|
|||
|
||||
super.sanityChecks(equations, t0, y0, t, y);
|
||||
|
||||
if ((vecAbsoluteTolerance != null) && (vecAbsoluteTolerance.length != y0.length)) {
|
||||
throw new IntegratorException(
|
||||
LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, y0.length, vecAbsoluteTolerance.length);
|
||||
if (equations instanceof ExtendedFirstOrderDifferentialEquations) {
|
||||
mainSetDimension = ((ExtendedFirstOrderDifferentialEquations) equations).getMainSetDimension();
|
||||
} else {
|
||||
mainSetDimension = equations.getDimension();
|
||||
}
|
||||
|
||||
if ((vecRelativeTolerance != null) && (vecRelativeTolerance.length != y0.length)) {
|
||||
if ((vecAbsoluteTolerance != null) && (vecAbsoluteTolerance.length != mainSetDimension)) {
|
||||
throw new IntegratorException(
|
||||
LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, y0.length, vecRelativeTolerance.length);
|
||||
LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, mainSetDimension, vecAbsoluteTolerance.length);
|
||||
}
|
||||
|
||||
if ((vecRelativeTolerance != null) && (vecRelativeTolerance.length != mainSetDimension)) {
|
||||
throw new IntegratorException(
|
||||
LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, mainSetDimension, vecRelativeTolerance.length);
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -187,7 +205,7 @@ public abstract class AdaptiveStepsizeIntegrator
|
|||
* @param equations differential equations set
|
||||
* @param forward forward integration indicator
|
||||
* @param order order of the method
|
||||
* @param scale scaling vector for the state vector
|
||||
* @param scale scaling vector for the state vector (can be shorter than state vector)
|
||||
* @param t0 start time
|
||||
* @param y0 state vector at t0
|
||||
* @param yDot0 first time derivative of y0
|
||||
|
|
@ -213,7 +231,7 @@ public abstract class AdaptiveStepsizeIntegrator
|
|||
double ratio;
|
||||
double yOnScale2 = 0;
|
||||
double yDotOnScale2 = 0;
|
||||
for (int j = 0; j < y0.length; ++j) {
|
||||
for (int j = 0; j < scale.length; ++j) {
|
||||
ratio = y0[j] / scale[j];
|
||||
yOnScale2 += ratio * ratio;
|
||||
ratio = yDot0[j] / scale[j];
|
||||
|
|
@ -234,7 +252,7 @@ public abstract class AdaptiveStepsizeIntegrator
|
|||
|
||||
// estimate the second derivative of the solution
|
||||
double yDDotOnScale = 0;
|
||||
for (int j = 0; j < y0.length; ++j) {
|
||||
for (int j = 0; j < scale.length; ++j) {
|
||||
ratio = (yDot1[j] - yDot0[j]) / scale[j];
|
||||
yDDotOnScale += ratio * ratio;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -135,7 +135,7 @@ public class DormandPrince54Integrator extends EmbeddedRungeKuttaIntegrator {
|
|||
|
||||
double error = 0;
|
||||
|
||||
for (int j = 0; j < y0.length; ++j) {
|
||||
for (int j = 0; j < mainSetDimension; ++j) {
|
||||
final double errSum = E1 * yDotK[0][j] + E3 * yDotK[2][j] +
|
||||
E4 * yDotK[3][j] + E5 * yDotK[4][j] +
|
||||
E6 * yDotK[5][j] + E7 * yDotK[6][j];
|
||||
|
|
@ -149,7 +149,7 @@ public class DormandPrince54Integrator extends EmbeddedRungeKuttaIntegrator {
|
|||
|
||||
}
|
||||
|
||||
return Math.sqrt(error / y0.length);
|
||||
return Math.sqrt(error / mainSetDimension);
|
||||
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -249,7 +249,7 @@ public class DormandPrince853Integrator extends EmbeddedRungeKuttaIntegrator {
|
|||
double error1 = 0;
|
||||
double error2 = 0;
|
||||
|
||||
for (int j = 0; j < y0.length; ++j) {
|
||||
for (int j = 0; j < mainSetDimension; ++j) {
|
||||
final double errSum1 = E1_01 * yDotK[0][j] + E1_06 * yDotK[5][j] +
|
||||
E1_07 * yDotK[6][j] + E1_08 * yDotK[7][j] +
|
||||
E1_09 * yDotK[8][j] + E1_10 * yDotK[9][j] +
|
||||
|
|
@ -274,7 +274,7 @@ public class DormandPrince853Integrator extends EmbeddedRungeKuttaIntegrator {
|
|||
den = 1.0;
|
||||
}
|
||||
|
||||
return Math.abs(h) * error1 / Math.sqrt(y0.length * den);
|
||||
return Math.abs(h) * error1 / Math.sqrt(mainSetDimension * den);
|
||||
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -242,7 +242,7 @@ public abstract class EmbeddedRungeKuttaIntegrator
|
|||
}
|
||||
|
||||
if (firstTime) {
|
||||
final double[] scale = new double[y0.length];
|
||||
final double[] scale = new double[mainSetDimension];
|
||||
if (vecAbsoluteTolerance == null) {
|
||||
for (int i = 0; i < scale.length; ++i) {
|
||||
scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * Math.abs(y[i]);
|
||||
|
|
|
|||
|
|
@ -429,7 +429,7 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
|
|||
/** Update scaling array.
|
||||
* @param y1 first state vector to use for scaling
|
||||
* @param y2 second state vector to use for scaling
|
||||
* @param scale scaling array to update
|
||||
* @param scale scaling array to update (can be shorter than state)
|
||||
*/
|
||||
private void rescale(final double[] y1, final double[] y2, final double[] scale) {
|
||||
if (vecAbsoluteTolerance == null) {
|
||||
|
|
@ -451,7 +451,7 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
|
|||
* @param y0 initial value of the state vector at t0
|
||||
* @param step global step
|
||||
* @param k iteration number (from 0 to sequence.length - 1)
|
||||
* @param scale scaling array
|
||||
* @param scale scaling array (can be shorter than state)
|
||||
* @param f placeholder where to put the state vector derivatives at each substep
|
||||
* (element 0 already contains initial derivative)
|
||||
* @param yMiddle placeholder where to put the state vector at the middle of the step
|
||||
|
|
@ -500,12 +500,12 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
|
|||
// stability check
|
||||
if (performTest && (j <= maxChecks) && (k < maxIter)) {
|
||||
double initialNorm = 0.0;
|
||||
for (int l = 0; l < y0.length; ++l) {
|
||||
for (int l = 0; l < scale.length; ++l) {
|
||||
final double ratio = f[0][l] / scale[l];
|
||||
initialNorm += ratio * ratio;
|
||||
}
|
||||
double deltaNorm = 0.0;
|
||||
for (int l = 0; l < y0.length; ++l) {
|
||||
for (int l = 0; l < scale.length; ++l) {
|
||||
final double ratio = (f[j+1][l] - f[0][l]) / scale[l];
|
||||
deltaNorm += ratio * ratio;
|
||||
}
|
||||
|
|
@ -607,7 +607,7 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
|
|||
}
|
||||
|
||||
// initial scaling
|
||||
final double[] scale = new double[y0.length];
|
||||
final double[] scale = new double[mainSetDimension];
|
||||
rescale(y, y, scale);
|
||||
|
||||
// initial order selection
|
||||
|
|
@ -709,11 +709,11 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
|
|||
|
||||
// estimate the error at the end of the step.
|
||||
error = 0;
|
||||
for (int j = 0; j < y0.length; ++j) {
|
||||
for (int j = 0; j < mainSetDimension; ++j) {
|
||||
final double e = Math.abs(y1[j] - y1Diag[0][j]) / scale[j];
|
||||
error += e * e;
|
||||
}
|
||||
error = Math.sqrt(error / y0.length);
|
||||
error = Math.sqrt(error / mainSetDimension);
|
||||
|
||||
if ((error > 1.0e15) || ((k > 1) && (error > maxError))) {
|
||||
// error is too big, we reduce the global step
|
||||
|
|
|
|||
|
|
@ -297,11 +297,11 @@ class GraggBulirschStoerStepInterpolator
|
|||
public double estimateError(final double[] scale) {
|
||||
double error = 0;
|
||||
if (currentDegree >= 5) {
|
||||
for (int i = 0; i < currentState.length; ++i) {
|
||||
for (int i = 0; i < scale.length; ++i) {
|
||||
final double e = polynoms[currentDegree][i] / scale[i];
|
||||
error += e * e;
|
||||
}
|
||||
error = Math.sqrt(error / currentState.length) * errfac[currentDegree-5];
|
||||
error = Math.sqrt(error / scale.length) * errfac[currentDegree - 5];
|
||||
}
|
||||
return error;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -108,7 +108,7 @@ public class HighamHall54Integrator extends EmbeddedRungeKuttaIntegrator {
|
|||
|
||||
double error = 0;
|
||||
|
||||
for (int j = 0; j < y0.length; ++j) {
|
||||
for (int j = 0; j < mainSetDimension; ++j) {
|
||||
double errSum = STATIC_E[0] * yDotK[0][j];
|
||||
for (int l = 1; l < STATIC_E.length; ++l) {
|
||||
errSum += STATIC_E[l] * yDotK[l][j];
|
||||
|
|
@ -123,7 +123,7 @@ public class HighamHall54Integrator extends EmbeddedRungeKuttaIntegrator {
|
|||
|
||||
}
|
||||
|
||||
return Math.sqrt(error / y0.length);
|
||||
return Math.sqrt(error / mainSetDimension);
|
||||
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -52,6 +52,12 @@ The <action> type attribute can be add,update,fix,remove.
|
|||
If the output is not quite correct, check for invisible trailing spaces!
|
||||
-->
|
||||
<release version="2.2" date="TBD" description="TBD">
|
||||
<action dev="luc" type="add" issue="MATH-388">
|
||||
Added a feature allowing error estimation to be computed only on a subset of
|
||||
Ordinary Differential Equations, considered as the main set, the remaining equations
|
||||
being considered only as an extension set that should not influence the ODE integration
|
||||
algorithm
|
||||
</action>
|
||||
<action dev="erans" type="fix" issue="MATH-382">
|
||||
Fixed bug in precondition check (method "setMicrosphereElements").
|
||||
</action>
|
||||
|
|
|
|||
|
|
@ -40,7 +40,7 @@ public class FirstOrderIntegratorWithJacobiansTest {
|
|||
// 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, 1.0e-4, 1.0e-4);
|
||||
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();
|
||||
|
|
@ -64,7 +64,7 @@ public class FirstOrderIntegratorWithJacobiansTest {
|
|||
public void testHighAccuracyExternalDifferentiation()
|
||||
throws IntegratorException, DerivativeException {
|
||||
FirstOrderIntegrator integ =
|
||||
new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10);
|
||||
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();
|
||||
|
|
@ -92,7 +92,7 @@ public class FirstOrderIntegratorWithJacobiansTest {
|
|||
public void testInternalDifferentiation()
|
||||
throws IntegratorException, DerivativeException {
|
||||
FirstOrderIntegrator integ =
|
||||
new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-4, 1.0e-4);
|
||||
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();
|
||||
|
|
@ -109,23 +109,23 @@ public class FirstOrderIntegratorWithJacobiansTest {
|
|||
extInt.setMaxEvaluations(5000);
|
||||
extInt.integrate(0, z, new double[][] { { 0.0 }, { 1.0 } }, 20.0, z, dZdZ0, dZdP);
|
||||
Assert.assertEquals(5000, extInt.getMaxEvaluations());
|
||||
Assert.assertTrue(extInt.getEvaluations() > 2000);
|
||||
Assert.assertTrue(extInt.getEvaluations() < 2500);
|
||||
Assert.assertTrue(extInt.getEvaluations() > 1500);
|
||||
Assert.assertTrue(extInt.getEvaluations() < 2100);
|
||||
Assert.assertEquals(4 * integ.getEvaluations(), extInt.getEvaluations());
|
||||
residualsP0.addValue(dZdP[0][0] - brusselator.dYdP0());
|
||||
residualsP1.addValue(dZdP[1][0] - brusselator.dYdP1());
|
||||
}
|
||||
Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.006);
|
||||
Assert.assertTrue(residualsP0.getStandardDeviation() < 0.0009);
|
||||
Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.009);
|
||||
Assert.assertTrue(residualsP1.getStandardDeviation() < 0.0014);
|
||||
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()
|
||||
throws IntegratorException, DerivativeException {
|
||||
FirstOrderIntegrator integ =
|
||||
new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-4, 1.0e-4);
|
||||
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) {
|
||||
|
|
@ -139,22 +139,22 @@ public class FirstOrderIntegratorWithJacobiansTest {
|
|||
extInt.setMaxEvaluations(5000);
|
||||
extInt.integrate(0, z, new double[][] { { 0.0 }, { 1.0 } }, 20.0, z, dZdZ0, dZdP);
|
||||
Assert.assertEquals(5000, extInt.getMaxEvaluations());
|
||||
Assert.assertTrue(extInt.getEvaluations() > 510);
|
||||
Assert.assertTrue(extInt.getEvaluations() < 610);
|
||||
Assert.assertTrue(extInt.getEvaluations() > 350);
|
||||
Assert.assertTrue(extInt.getEvaluations() < 510);
|
||||
Assert.assertEquals(integ.getEvaluations(), extInt.getEvaluations());
|
||||
residualsP0.addValue(dZdP[0][0] - brusselator.dYdP0());
|
||||
residualsP1.addValue(dZdP[1][0] - brusselator.dYdP1());
|
||||
}
|
||||
Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.004);
|
||||
Assert.assertTrue(residualsP0.getStandardDeviation() < 0.0008);
|
||||
Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.005);
|
||||
Assert.assertTrue(residualsP1.getStandardDeviation() < 0.0010);
|
||||
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() throws IntegratorException, DerivativeException {
|
||||
FirstOrderIntegrator integ =
|
||||
new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10);
|
||||
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);
|
||||
double[][] dydy0 = new double[2][2];
|
||||
|
|
@ -164,16 +164,16 @@ public class FirstOrderIntegratorWithJacobiansTest {
|
|||
new FirstOrderIntegratorWithJacobians(integ, circle);
|
||||
extInt.integrate(0, y, circle.exactDyDp(0), t, y, dydy0, dydp);
|
||||
for (int i = 0; i < y.length; ++i) {
|
||||
Assert.assertEquals(circle.exactY(t)[i], y[i], 1.0e-10);
|
||||
Assert.assertEquals(circle.exactY(t)[i], y[i], 1.0e-9);
|
||||
}
|
||||
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-10);
|
||||
Assert.assertEquals(circle.exactDyDy0(t)[i][j], dydy0[i][j], 1.0e-9);
|
||||
}
|
||||
}
|
||||
for (int i = 0; i < dydp.length; ++i) {
|
||||
for (int j = 0; j < dydp[i].length; ++j) {
|
||||
Assert.assertEquals(circle.exactDyDp(t)[i][j], dydp[i][j], 1.0e-8);
|
||||
Assert.assertEquals(circle.exactDyDp(t)[i][j], dydp[i][j], 1.0e-7);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -181,7 +181,7 @@ public class FirstOrderIntegratorWithJacobiansTest {
|
|||
@Test
|
||||
public void testStepHandlerResult() throws IntegratorException, DerivativeException {
|
||||
FirstOrderIntegrator integ =
|
||||
new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10);
|
||||
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 };
|
||||
final Circle circle = new Circle(y, 1.0, 1.0, 0.1);
|
||||
double[][] dydy0 = new double[2][2];
|
||||
|
|
@ -207,16 +207,16 @@ public class FirstOrderIntegratorWithJacobiansTest {
|
|||
Assert.assertEquals(interpolator.getPreviousTime(), extInt.getCurrentStepStart(), 1.0e-10);
|
||||
Assert.assertTrue(extInt.getCurrentSignedStepsize() < 0.5);
|
||||
for (int i = 0; i < y.length; ++i) {
|
||||
Assert.assertEquals(circle.exactY(t)[i], y[i], 1.0e-10);
|
||||
Assert.assertEquals(circle.exactY(t)[i], y[i], 1.0e-9);
|
||||
}
|
||||
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-10);
|
||||
Assert.assertEquals(circle.exactDyDy0(t)[i][j], dydy0[i][j], 1.0e-9);
|
||||
}
|
||||
}
|
||||
for (int i = 0; i < dydp.length; ++i) {
|
||||
for (int j = 0; j < dydp[i].length; ++j) {
|
||||
Assert.assertEquals(circle.exactDyDp(t)[i][j], dydp[i][j], 1.0e-8);
|
||||
Assert.assertEquals(circle.exactDyDp(t)[i][j], dydp[i][j], 3.0e-8);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -225,16 +225,16 @@ public class FirstOrderIntegratorWithJacobiansTest {
|
|||
double[][] dydpDot = interpolator.getInterpolatedDyDpDot();
|
||||
|
||||
for (int i = 0; i < yDot.length; ++i) {
|
||||
Assert.assertEquals(circle.exactYDot(t)[i], yDot[i], 1.0e-11);
|
||||
Assert.assertEquals(circle.exactYDot(t)[i], yDot[i], 1.0e-10);
|
||||
}
|
||||
for (int i = 0; i < dydy0Dot.length; ++i) {
|
||||
for (int j = 0; j < dydy0Dot[i].length; ++j) {
|
||||
Assert.assertEquals(circle.exactDyDy0Dot(t)[i][j], dydy0Dot[i][j], 1.0e-11);
|
||||
Assert.assertEquals(circle.exactDyDy0Dot(t)[i][j], dydy0Dot[i][j], 1.0e-10);
|
||||
}
|
||||
}
|
||||
for (int i = 0; i < dydpDot.length; ++i) {
|
||||
for (int j = 0; j < dydpDot[i].length; ++j) {
|
||||
Assert.assertEquals(circle.exactDyDpDot(t)[i][j], dydpDot[i][j], 1.0e-9);
|
||||
Assert.assertEquals(circle.exactDyDpDot(t)[i][j], dydpDot[i][j], 3.0e-9);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -245,7 +245,7 @@ public class FirstOrderIntegratorWithJacobiansTest {
|
|||
@Test
|
||||
public void testEventHandler() throws IntegratorException, DerivativeException {
|
||||
FirstOrderIntegrator integ =
|
||||
new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10);
|
||||
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 };
|
||||
final Circle circle = new Circle(y, 1.0, 1.0, 0.1);
|
||||
double[][] dydy0 = new double[2][2];
|
||||
|
|
|
|||
Loading…
Reference in New Issue