diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/ode/nonstiff/AdaptiveStepsizeFieldIntegrator.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math3/ode/nonstiff/AdaptiveStepsizeFieldIntegrator.java | 366 |
1 files changed, 366 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/AdaptiveStepsizeFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdaptiveStepsizeFieldIntegrator.java new file mode 100644 index 0000000..c8e592b --- /dev/null +++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdaptiveStepsizeFieldIntegrator.java @@ -0,0 +1,366 @@ +/* + * 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.math3.ode.nonstiff; + +import org.apache.commons.math3.Field; +import org.apache.commons.math3.RealFieldElement; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.ode.AbstractFieldIntegrator; +import org.apache.commons.math3.ode.FieldEquationsMapper; +import org.apache.commons.math3.ode.FieldODEState; +import org.apache.commons.math3.ode.FieldODEStateAndDerivative; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.MathArrays; +import org.apache.commons.math3.util.MathUtils; + +/** + * This abstract class holds the common part of all adaptive + * stepsize integrators for Ordinary Differential Equations. + * + * <p>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 + * <pre> + * threshold_i = absTol_i + relTol_i * max (abs (ym), abs (ym+1)) + * </pre> + * 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> + * <p> + * Note that <em>only</em> the {@link FieldODEState#getState() main part} + * of the state vector is used for stepsize control. The {@link + * FieldODEState#getSecondaryState(int) secondary parts} of the state + * vector are explicitly ignored for stepsize control. + * </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 main set dimension) then the step is accepted, + * otherwise the step is rejected and a new attempt is made with a new + * stepsize.</p> + * + * @param <T> the type of the field elements + * @since 3.6 + * + */ + +public abstract class AdaptiveStepsizeFieldIntegrator<T extends RealFieldElement<T>> + extends AbstractFieldIntegrator<T> { + + /** 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; + + /** Main set dimension. */ + protected int mainSetDimension; + + /** User supplied initial step. */ + private T initialStep; + + /** Minimal step. */ + private T minStep; + + /** Maximal step. */ + private T maxStep; + + /** Build an integrator with the given stepsize bounds. + * The default step handler does nothing. + * @param field field to which the time and state vector elements belong + * @param name name of the method + * @param minStep minimal step (sign is irrelevant, regardless of + * integration direction, forward or backward), the last step can + * be smaller than this + * @param maxStep maximal step (sign is irrelevant, regardless of + * integration direction, forward or backward), the last step can + * be smaller than this + * @param scalAbsoluteTolerance allowed absolute error + * @param scalRelativeTolerance allowed relative error + */ + public AdaptiveStepsizeFieldIntegrator(final Field<T> field, final String name, + final double minStep, final double maxStep, + final double scalAbsoluteTolerance, + final double scalRelativeTolerance) { + + super(field, name); + setStepSizeControl(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); + resetInternalState(); + + } + + /** Build an integrator with the given stepsize bounds. + * The default step handler does nothing. + * @param field field to which the time and state vector elements belong + * @param name name of the method + * @param minStep minimal step (sign is irrelevant, regardless of + * integration direction, forward or backward), the last step can + * be smaller than this + * @param maxStep maximal step (sign is irrelevant, regardless of + * integration direction, forward or backward), the last step can + * be smaller than this + * @param vecAbsoluteTolerance allowed absolute error + * @param vecRelativeTolerance allowed relative error + */ + public AdaptiveStepsizeFieldIntegrator(final Field<T> field, final String name, + final double minStep, final double maxStep, + final double[] vecAbsoluteTolerance, + final double[] vecRelativeTolerance) { + + super(field, name); + setStepSizeControl(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); + resetInternalState(); + + } + + /** Set the adaptive step size control parameters. + * <p> + * A side effect of this method is to also reset the initial + * step so it will be automatically computed by the integrator + * if {@link #setInitialStepSize(RealFieldElement) setInitialStepSize} + * is not called by the user. + * </p> + * @param minimalStep minimal step (must be positive even for backward + * integration), the last step can be smaller than this + * @param maximalStep maximal step (must be positive even for backward + * integration) + * @param absoluteTolerance allowed absolute error + * @param relativeTolerance allowed relative error + */ + public void setStepSizeControl(final double minimalStep, final double maximalStep, + final double absoluteTolerance, + final double relativeTolerance) { + + minStep = getField().getZero().add(FastMath.abs(minimalStep)); + maxStep = getField().getZero().add(FastMath.abs(maximalStep)); + initialStep = getField().getOne().negate(); + + scalAbsoluteTolerance = absoluteTolerance; + scalRelativeTolerance = relativeTolerance; + vecAbsoluteTolerance = null; + vecRelativeTolerance = null; + + } + + /** Set the adaptive step size control parameters. + * <p> + * A side effect of this method is to also reset the initial + * step so it will be automatically computed by the integrator + * if {@link #setInitialStepSize(RealFieldElement) setInitialStepSize} + * is not called by the user. + * </p> + * @param minimalStep minimal step (must be positive even for backward + * integration), the last step can be smaller than this + * @param maximalStep maximal step (must be positive even for backward + * integration) + * @param absoluteTolerance allowed absolute error + * @param relativeTolerance allowed relative error + */ + public void setStepSizeControl(final double minimalStep, final double maximalStep, + final double[] absoluteTolerance, + final double[] relativeTolerance) { + + minStep = getField().getZero().add(FastMath.abs(minimalStep)); + maxStep = getField().getZero().add(FastMath.abs(maximalStep)); + initialStep = getField().getOne().negate(); + + scalAbsoluteTolerance = 0; + scalRelativeTolerance = 0; + vecAbsoluteTolerance = absoluteTolerance.clone(); + vecRelativeTolerance = relativeTolerance.clone(); + + } + + /** Set the initial step size. + * <p>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.</p> + * @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 T initialStepSize) { + if (initialStepSize.subtract(minStep).getReal() < 0 || + initialStepSize.subtract(maxStep).getReal() > 0) { + initialStep = getField().getOne().negate(); + } else { + initialStep = initialStepSize; + } + } + + /** {@inheritDoc} */ + @Override + protected void sanityChecks(final FieldODEState<T> eqn, final T t) + throws DimensionMismatchException, NumberIsTooSmallException { + + super.sanityChecks(eqn, t); + + mainSetDimension = eqn.getStateDimension(); + + if (vecAbsoluteTolerance != null && vecAbsoluteTolerance.length != mainSetDimension) { + throw new DimensionMismatchException(mainSetDimension, vecAbsoluteTolerance.length); + } + + if (vecRelativeTolerance != null && vecRelativeTolerance.length != mainSetDimension) { + throw new DimensionMismatchException(mainSetDimension, vecRelativeTolerance.length); + } + + } + + /** Initialize the integration step. + * @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) + * @param state0 state at integration start time + * @param mapper mapper for all the equations + * @return first integration step + * @exception MaxCountExceededException if the number of functions evaluations is exceeded + * @exception DimensionMismatchException if arrays dimensions do not match equations settings + */ + public T initializeStep(final boolean forward, final int order, final T[] scale, + final FieldODEStateAndDerivative<T> state0, + final FieldEquationsMapper<T> mapper) + throws MaxCountExceededException, DimensionMismatchException { + + if (initialStep.getReal() > 0) { + // use the user provided value + return forward ? initialStep : initialStep.negate(); + } + + // very rough first guess : h = 0.01 * ||y/scale|| / ||y'/scale|| + // this guess will be used to perform an Euler step + final T[] y0 = mapper.mapState(state0); + final T[] yDot0 = mapper.mapDerivative(state0); + T yOnScale2 = getField().getZero(); + T yDotOnScale2 = getField().getZero(); + for (int j = 0; j < scale.length; ++j) { + final T ratio = y0[j].divide(scale[j]); + yOnScale2 = yOnScale2.add(ratio.multiply(ratio)); + final T ratioDot = yDot0[j].divide(scale[j]); + yDotOnScale2 = yDotOnScale2.add(ratioDot.multiply(ratioDot)); + } + + T h = (yOnScale2.getReal() < 1.0e-10 || yDotOnScale2.getReal() < 1.0e-10) ? + getField().getZero().add(1.0e-6) : + yOnScale2.divide(yDotOnScale2).sqrt().multiply(0.01); + if (! forward) { + h = h.negate(); + } + + // perform an Euler step using the preceding rough guess + final T[] y1 = MathArrays.buildArray(getField(), y0.length); + for (int j = 0; j < y0.length; ++j) { + y1[j] = y0[j].add(yDot0[j].multiply(h)); + } + final T[] yDot1 = computeDerivatives(state0.getTime().add(h), y1); + + // estimate the second derivative of the solution + T yDDotOnScale = getField().getZero(); + for (int j = 0; j < scale.length; ++j) { + final T ratioDotDot = yDot1[j].subtract(yDot0[j]).divide(scale[j]); + yDDotOnScale = yDDotOnScale.add(ratioDotDot.multiply(ratioDotDot)); + } + yDDotOnScale = yDDotOnScale.sqrt().divide(h); + + // step size is computed such that + // h^order * max (||y'/tol||, ||y''/tol||) = 0.01 + final T maxInv2 = MathUtils.max(yDotOnScale2.sqrt(), yDDotOnScale); + final T h1 = maxInv2.getReal() < 1.0e-15 ? + MathUtils.max(getField().getZero().add(1.0e-6), h.abs().multiply(0.001)) : + maxInv2.multiply(100).reciprocal().pow(1.0 / order); + h = MathUtils.min(h.abs().multiply(100), h1); + h = MathUtils.max(h, state0.getTime().abs().multiply(1.0e-12)); // avoids cancellation when computing t1 - t0 + h = MathUtils.max(minStep, MathUtils.min(maxStep, h)); + if (! forward) { + h = h.negate(); + } + + 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 NumberIsTooSmallException if the step is too small and acceptSmall is false + */ + protected T filterStep(final T h, final boolean forward, final boolean acceptSmall) + throws NumberIsTooSmallException { + + T filteredH = h; + if (h.abs().subtract(minStep).getReal() < 0) { + if (acceptSmall) { + filteredH = forward ? minStep : minStep.negate(); + } else { + throw new NumberIsTooSmallException(LocalizedFormats.MINIMAL_STEPSIZE_REACHED_DURING_INTEGRATION, + h.abs().getReal(), minStep.getReal(), true); + } + } + + if (filteredH.subtract(maxStep).getReal() > 0) { + filteredH = maxStep; + } else if (filteredH.add(maxStep).getReal() < 0) { + filteredH = maxStep.negate(); + } + + return filteredH; + + } + + /** Reset internal state to dummy values. */ + protected void resetInternalState() { + setStepStart(null); + setStepSize(minStep.multiply(maxStep).sqrt()); + } + + /** Get the minimal step. + * @return minimal step + */ + public T getMinStep() { + return minStep; + } + + /** Get the maximal step. + * @return maximal step + */ + public T getMaxStep() { + return maxStep; + } + +} |