org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of commons-math Show documentation
Show all versions of commons-math Show documentation
The Math project is a library of lightweight, self-contained mathematics and statistics components addressing the most common practical problems not immediately available in the Java programming language or commons-lang.
/*
* 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.nonstiff;
import org.apache.commons.math.ode.AbstractIntegrator;
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.IntegratorException;
/**
* This abstract class holds the common part of all adaptive
* stepsize integrators for Ordinary Differential Equations.
*
* These algorithms perform integration with stepsize control, which
* means the user does not specify the integration step but rather a
* tolerance on error. The error threshold is computed as
*
* threshold_i = absTol_i + relTol_i * max (abs (ym), abs (ym+1))
*
* 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.
*
* If the estimated error for ym+1 is such that
*
* sqrt((sum (errEst_i / threshold_i)^2 ) / n) < 1
*
*
* (where n is the state vector dimension) then the step is accepted,
* otherwise the step is rejected and a new attempt is made with a new
* stepsize.
*
* @version $Revision: 795591 $ $Date: 2009-07-19 14:36:46 -0400 (Sun, 19 Jul 2009) $
* @since 1.2
*
*/
public abstract class AdaptiveStepsizeIntegrator
extends AbstractIntegrator {
/** Build an integrator with the given stepsize bounds.
* The default step handler does nothing.
* @param name name of the method
* @param minStep minimal step (must be positive even for backward
* integration), the last step can be smaller than this
* @param maxStep maximal step (must be positive even for backward
* integration)
* @param scalAbsoluteTolerance allowed absolute error
* @param scalRelativeTolerance allowed relative error
*/
public AdaptiveStepsizeIntegrator(final String name,
final double minStep, final double maxStep,
final double scalAbsoluteTolerance,
final double scalRelativeTolerance) {
super(name);
this.minStep = Math.abs(minStep);
this.maxStep = Math.abs(maxStep);
this.initialStep = -1.0;
this.scalAbsoluteTolerance = scalAbsoluteTolerance;
this.scalRelativeTolerance = scalRelativeTolerance;
this.vecAbsoluteTolerance = null;
this.vecRelativeTolerance = null;
resetInternalState();
}
/** Build an integrator with the given stepsize bounds.
* The default step handler does nothing.
* @param name name of the method
* @param minStep minimal step (must be positive even for backward
* integration), the last step can be smaller than this
* @param maxStep maximal step (must be positive even for backward
* integration)
* @param vecAbsoluteTolerance allowed absolute error
* @param vecRelativeTolerance allowed relative error
*/
public AdaptiveStepsizeIntegrator(final String name,
final double minStep, final double maxStep,
final double[] vecAbsoluteTolerance,
final double[] vecRelativeTolerance) {
super(name);
this.minStep = minStep;
this.maxStep = maxStep;
this.initialStep = -1.0;
this.scalAbsoluteTolerance = 0;
this.scalRelativeTolerance = 0;
this.vecAbsoluteTolerance = vecAbsoluteTolerance.clone();
this.vecRelativeTolerance = vecRelativeTolerance.clone();
resetInternalState();
}
/** Set the initial step size.
* This method allows the user to specify an initial positive
* step size instead of letting the integrator guess it by
* itself. If this method is not called before integration is
* started, the initial step size will be estimated by the
* integrator.
* @param initialStepSize initial step size to use (must be positive even
* for backward integration ; providing a negative value or a value
* outside of the min/max step interval will lead the integrator to
* ignore the value and compute the initial step size by itself)
*/
public void setInitialStepSize(final double initialStepSize) {
if ((initialStepSize < minStep) || (initialStepSize > maxStep)) {
initialStep = -1.0;
} else {
initialStep = initialStepSize;
}
}
/** Perform some sanity checks on the integration parameters.
* @param equations differential equations set
* @param t0 start time
* @param y0 state vector at t0
* @param t target time for the integration
* @param y placeholder where to put the state vector
* @exception IntegratorException if some inconsistency is detected
*/
@Override
protected void sanityChecks(final FirstOrderDifferentialEquations equations,
final double t0, final double[] y0,
final double t, final double[] y)
throws IntegratorException {
super.sanityChecks(equations, t0, y0, t, y);
if ((vecAbsoluteTolerance != null) && (vecAbsoluteTolerance.length != y0.length)) {
throw new IntegratorException(
"dimensions mismatch: state vector has dimension {0}," +
" absolute tolerance vector has dimension {1}",
y0.length, vecAbsoluteTolerance.length);
}
if ((vecRelativeTolerance != null) && (vecRelativeTolerance.length != y0.length)) {
throw new IntegratorException(
"dimensions mismatch: state vector has dimension {0}," +
" relative tolerance vector has dimension {1}",
y0.length, vecRelativeTolerance.length);
}
}
/** 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
* @param t0 start time
* @param y0 state vector at t0
* @param yDot0 first time derivative of y0
* @param y1 work array for a state vector
* @param yDot1 work array for the first time derivative of y1
* @return first integration step
* @exception DerivativeException this exception is propagated to
* the caller if the underlying user function triggers one
*/
public double initializeStep(final FirstOrderDifferentialEquations equations,
final boolean forward, final int order, final double[] scale,
final double t0, final double[] y0, final double[] yDot0,
final double[] y1, final double[] yDot1)
throws DerivativeException {
if (initialStep > 0) {
// use the user provided value
return forward ? initialStep : -initialStep;
}
// very rough first guess : h = 0.01 * ||y/scale|| / ||y'/scale||
// this guess will be used to perform an Euler step
double ratio;
double yOnScale2 = 0;
double yDotOnScale2 = 0;
for (int j = 0; j < y0.length; ++j) {
ratio = y0[j] / scale[j];
yOnScale2 += ratio * ratio;
ratio = yDot0[j] / scale[j];
yDotOnScale2 += ratio * ratio;
}
double h = ((yOnScale2 < 1.0e-10) || (yDotOnScale2 < 1.0e-10)) ?
1.0e-6 : (0.01 * Math.sqrt(yOnScale2 / yDotOnScale2));
if (! forward) {
h = -h;
}
// perform an Euler step using the preceding rough guess
for (int j = 0; j < y0.length; ++j) {
y1[j] = y0[j] + h * yDot0[j];
}
computeDerivatives(t0 + h, y1, yDot1);
// estimate the second derivative of the solution
double yDDotOnScale = 0;
for (int j = 0; j < y0.length; ++j) {
ratio = (yDot1[j] - yDot0[j]) / scale[j];
yDDotOnScale += ratio * ratio;
}
yDDotOnScale = Math.sqrt(yDDotOnScale) / h;
// step size is computed such that
// h^order * max (||y'/tol||, ||y''/tol||) = 0.01
final double maxInv2 = Math.max(Math.sqrt(yDotOnScale2), yDDotOnScale);
final double h1 = (maxInv2 < 1.0e-15) ?
Math.max(1.0e-6, 0.001 * Math.abs(h)) :
Math.pow(0.01 / maxInv2, 1.0 / order);
h = Math.min(100.0 * Math.abs(h), h1);
h = Math.max(h, 1.0e-12 * Math.abs(t0)); // avoids cancellation when computing t1 - t0
if (h < getMinStep()) {
h = getMinStep();
}
if (h > getMaxStep()) {
h = getMaxStep();
}
if (! forward) {
h = -h;
}
return h;
}
/** Filter the integration step.
* @param h signed step
* @param forward forward integration indicator
* @param acceptSmall if true, steps smaller than the minimal value
* are silently increased up to this value, if false such small
* steps generate an exception
* @return a bounded integration step (h if no bound is reach, or a bounded value)
* @exception IntegratorException if the step is too small and acceptSmall is false
*/
protected double filterStep(final double h, final boolean forward, final boolean acceptSmall)
throws IntegratorException {
double filteredH = h;
if (Math.abs(h) < minStep) {
if (acceptSmall) {
filteredH = forward ? minStep : -minStep;
} else {
throw new IntegratorException(
"minimal step size ({0}) reached, integration needs {1}",
minStep, Math.abs(h));
}
}
if (filteredH > maxStep) {
filteredH = maxStep;
} else if (filteredH < -maxStep) {
filteredH = -maxStep;
}
return filteredH;
}
/** {@inheritDoc} */
public abstract double integrate (FirstOrderDifferentialEquations equations,
double t0, double[] y0,
double t, double[] y)
throws DerivativeException, IntegratorException;
/** {@inheritDoc} */
@Override
public double getCurrentStepStart() {
return stepStart;
}
/** Reset internal state to dummy values. */
protected void resetInternalState() {
stepStart = Double.NaN;
stepSize = Math.sqrt(minStep * maxStep);
}
/** Get the minimal step.
* @return minimal step
*/
public double getMinStep() {
return minStep;
}
/** Get the maximal step.
* @return maximal step
*/
public double getMaxStep() {
return maxStep;
}
/** Minimal step. */
private double minStep;
/** Maximal step. */
private double maxStep;
/** User supplied initial step. */
private double initialStep;
/** Allowed absolute scalar error. */
protected double scalAbsoluteTolerance;
/** Allowed relative scalar error. */
protected double scalRelativeTolerance;
/** Allowed absolute vectorial error. */
protected double[] vecAbsoluteTolerance;
/** Allowed relative vectorial error. */
protected double[] vecRelativeTolerance;
}