adaptivestepsizeintegrator.java

Apache的common math数学软件包

/*
 * 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 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>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,
 * otherwise the step is rejected and a new attempt is made with a new
 * stepsize.</p>
 *
 * @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $
 * @since 1.2
 *
 */

public abstract class AdaptiveStepsizeIntegrator
  implements FirstOrderIntegrator {

  /** Build an integrator with the given stepsize bounds.
   * The default step handler does nothing.
   * @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(double minStep, double maxStep,
                                    double scalAbsoluteTolerance,
                                    double scalRelativeTolerance) {

    this.minStep     = minStep;
    this.maxStep     = maxStep;
    this.initialStep = -1.0;

    this.scalAbsoluteTolerance = scalAbsoluteTolerance;
    this.scalRelativeTolerance = scalRelativeTolerance;
    this.vecAbsoluteTolerance  = null;
    this.vecRelativeTolerance  = null;

    // set the default step handler
    handler = DummyStepHandler.getInstance();

    switchesHandler = new SwitchingFunctionsHandler();

    resetInternalState();

  }

  /** Build an integrator with the given stepsize bounds.
   * The default step handler does nothing.
   * @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(double minStep, double maxStep,
                                    double[] vecAbsoluteTolerance,
                                    double[] vecRelativeTolerance) {

    this.minStep     = minStep;
    this.maxStep     = maxStep;
    this.initialStep = -1.0;

    this.scalAbsoluteTolerance = 0;
    this.scalRelativeTolerance = 0;
    this.vecAbsoluteTolerance  = vecAbsoluteTolerance;
    this.vecRelativeTolerance  = vecRelativeTolerance;

    // set the default step handler
    handler = DummyStepHandler.getInstance();

    switchesHandler = new SwitchingFunctionsHandler();

    resetInternalState();

  }

  /** 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(double initialStepSize) {
    if ((initialStepSize < minStep) || (initialStepSize > maxStep)) {
      initialStep = -1.0;
    } else {
      initialStep = initialStepSize;
    }
  }

  /** Set the step handler for this integrator.
   * The handler will be called by the integrator for each accepted
   * step.
   * @param handler handler for the accepted steps
   */
  public void setStepHandler (StepHandler handler) {
    this.handler = handler;
  }

  /** Get the step handler for this integrator.
   * @return the step handler for this integrator
   */
  public StepHandler getStepHandler() {
    return handler;
  }

  /** Add a switching function to the integrator.
   * @param function switching function
   * @param maxCheckInterval maximal time interval between switching
   * function checks (this interval prevents missing sign changes in
   * case the integration steps becomes very large)
   * @param convergence convergence threshold in the event time search
   * @param maxIterationCount upper limit of the iteration count in
   * the event time search
   */
  public void addSwitchingFunction(SwitchingFunction function,
                                   double maxCheckInterval,
                                   double convergence,
                                   int maxIterationCount) {
    switchesHandler.add(function, maxCheckInterval, convergence, maxIterationCount);
  }

  /** 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
   */
  protected void sanityChecks(FirstOrderDifferentialEquations equations,
                              double t0, double[] y0, double t, double[] y)
    throws IntegratorException {
      if (equations.getDimension() != y0.length) {
          throw new IntegratorException("dimensions mismatch: ODE problem has dimension {0}," +
                                        " initial state vector has dimension {1}",
                                        new Object[] {
                                          new Integer(equations.getDimension()),
                                          new Integer(y0.length)
                                        });
      }
      if (equations.getDimension() != y.length) {
          throw new IntegratorException("dimensions mismatch: ODE problem has dimension {0}," +
                                        " final state vector has dimension {1}",
                                        new Object[] {
                                          new Integer(equations.getDimension()),
                                          new Integer(y.length)
                                        });
      }
      if ((vecAbsoluteTolerance != null) && (vecAbsoluteTolerance.length != y0.length)) {
          throw new IntegratorException("dimensions mismatch: state vector has dimension {0}," +
                                        " absolute tolerance vector has dimension {1}",
                                        new Object[] {
                                          new Integer(y0.length),
                                          new Integer(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}",
                                        new Object[] {
                                          new Integer(y0.length),
                                          new Integer(vecRelativeTolerance.length)
                                        });
      }
      if (Math.abs(t - t0) <= 1.0e-12 * Math.max(Math.abs(t0), Math.abs(t))) {
        throw new IntegratorException("too small integration interval: length = {0}",
                                      new Object[] { new Double(Math.abs(t - t0)) });
      }
      
  }

  /** Initialize the integration step.