dormandprince853stepinterpolator.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;

import java.io.ObjectOutput;
import java.io.ObjectInput;
import java.io.IOException;

/**
 * This class represents an interpolator over the last step during an
 * ODE integration for the 8(5,3) Dormand-Prince integrator.
 *
 * @see DormandPrince853Integrator
 *
 * @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $
 * @since 1.2
 */

class DormandPrince853StepInterpolator
  extends RungeKuttaStepInterpolator {

  /** Simple constructor.
   * This constructor builds an instance that is not usable yet, the
   * {@link #reinitialize} method should be called before using the
   * instance in order to initialize the internal arrays. This
   * constructor is used only in order to delay the initialization in
   * some cases. The {@link EmbeddedRungeKuttaIntegrator} uses the
   * prototyping design pattern to create the step interpolators by
   * cloning an uninitialized model and latter initializing the copy.
   */
  public DormandPrince853StepInterpolator() {
    super();
    yDotKLast = null;
    v         = null;
    vectorsInitialized = false;
  }

  /** Copy constructor.
   * @param interpolator interpolator to copy from. The copy is a deep
   * copy: its arrays are separated from the original arrays of the
   * instance
   */
  public DormandPrince853StepInterpolator(DormandPrince853StepInterpolator interpolator) {

    super(interpolator);

    if (interpolator.currentState == null) {

      yDotKLast = null;
      v         = null;
      vectorsInitialized = false;

    } else {

      int dimension = interpolator.currentState.length;

      yDotKLast    = new double[3][];
      for (int k = 0; k < yDotKLast.length; ++k) {
        yDotKLast[k] = new double[dimension];
        System.arraycopy(interpolator.yDotKLast[k], 0, yDotKLast[k], 0,
                         dimension);
      }

      v = new double[7][];
      for (int k = 0; k < v.length; ++k) {
        v[k] = new double[dimension];
        System.arraycopy(interpolator.v[k], 0, v[k], 0, dimension);
      }

      vectorsInitialized = interpolator.vectorsInitialized;

    }

  }

  /** Really copy the finalized instance.
   * @return a copy of the finalized instance
   */
  protected StepInterpolator doCopy() {
    return new DormandPrince853StepInterpolator(this);
  }

  /** Reinitialize the instance
   * Some embedded Runge-Kutta integrators need fewer functions
   * evaluations than their counterpart step interpolators. So the
   * interpolator should perform the last evaluations they need by
   * themselves. The {@link EmbeddedRungeKuttaIntegrator
   * EmbeddedRungeKuttaIntegrator} abstract class calls this method in
   * order to let the step interpolator perform the evaluations it
   * needs. These evaluations will be performed during the call to
   * <code>doFinalize</code> if any, i.e. only if the step handler
   * either calls the {@link AbstractStepInterpolator#finalizeStep
   * finalizeStep} method or the {@link
   * AbstractStepInterpolator#getInterpolatedState getInterpolatedState}
   * method (for an interpolator which needs a finalization) or if it clones
   * the step interpolator.
   * @param equations set of differential equations being integrated
   * @param y reference to the integrator array holding the state at
   * the end of the step
   * @param yDotK reference to the integrator array holding all the
   * intermediate slopes
   * @param forward integration direction indicator
   */
  public void reinitialize(FirstOrderDifferentialEquations equations,
                           double[] y, double[][] yDotK, boolean forward) {

    super.reinitialize(equations, y, yDotK, forward);

    int dimension = currentState.length;

    yDotKLast = new double[3][];
    for (int k = 0; k < yDotKLast.length; ++k) {
      yDotKLast[k] = new double[dimension];
    }

    v = new double[7][];
    for (int k = 0; k < v.length; ++k) {
      v[k]  = new double[dimension];
    }

    vectorsInitialized = false;

  }

  /** Store the current step time.
   * @param t current time
   */
  public void storeTime(double t) {
    super.storeTime(t);
    vectorsInitialized = false;
  }

  /** Compute the state at the interpolated time.
   * This is the main processing method that should be implemented by
   * the derived classes to perform the interpolation.
   * @param theta normalized interpolation abscissa within the step
   * (theta is zero at the previous time step and one at the current time step)
   * @param oneMinusThetaH time gap between the interpolated time and
   * the current time
   * @throws DerivativeException this exception is propagated to the caller if the
   * underlying user function triggers one
   */
  protected void computeInterpolatedState(double theta,
                                          double oneMinusThetaH)
    throws DerivativeException {

    if (! vectorsInitialized) {

      if (v == null) {
        v = new double[7][];
        for (int k = 0; k < 7; ++k) {
          v[k] = new double[interpolatedState.length];
        }
      }

      // perform the last evaluations if they have not been done yet
      finalizeStep();

      // compute the interpolation vectors for this time step
      for (int i = 0; i < interpolatedState.length; ++i) {
        v[0][i] = h * (b_01 * yDotK[0][i]  + b_06 * yDotK[5][i] + b_07 * yDotK[6][i] +
                       b_08 * yDotK[7][i]  + b_09 * yDotK[8][i] + b_10 * yDotK[9][i] +
                       b_11 * yDotK[10][i] + b_12 * yDotK[11][i]);
        v[1][i] = h * yDotK[0][i] - v[0][i];
        v[2][i] = v[0][i] - v[1][i] - h * yDotK[12][i];
        for (int k = 0; k < d.length; ++k) {
          v[k+3][i] = h * (d[k][0] * yDotK[0][i]  + d[k][1] * yDotK[5][i]  + d[k][2] * yDotK[6][i] +
                           d[k][3] * yDotK[7][i]  + d[k][4] * yDotK[8][i]  + d[k][5] * yDotK[9][i] +
                           d[k][6] * yDotK[10][i] + d[k][7] * yDotK[11][i] + d[k][8] * yDotK[12][i] +
                           d[k][9]  * yDotKLast[0][i] +
                           d[k][10] * yDotKLast[1][i] +
                           d[k][11] * yDotKLast[2][i]);
        }
      }

      vectorsInitialized = true;

    }

    double eta = oneMinusThetaH / h;

    for (int i = 0; i < interpolatedState.length; ++i) {
      interpolatedState[i] =
          currentState[i] - eta * (v[0][i] - theta * (v[1][i] +
                  theta * (v[2][i] + eta * (v[3][i] + theta * (v[4][i] +
                          eta * (v[5][i] + theta * (v[6][i])))))));
    }

  }
 
  /**
   * Really finalize the step.
   * Perform the last 3 functions evaluations (k14, k15, k16)
   * @throws DerivativeException this exception is propagated to the caller if the
   * underlying user function triggers one
   */
  protected void doFinalize()
    throws DerivativeException {

    if (currentState == null) {
      // we are finalizing an uninitialized instance
      return;
    }

    double s;
    double[] yTmp = new double[currentState.length];

    // k14
    for (int j = 0; j < currentState.length; ++j) {
      s = k14_01 * yDotK[0][j]  + k14_06 * yDotK[5][j]  + k14_07 * yDotK[6][j] +
          k14_08 * yDotK[7][j]  + k14_09 * yDotK[8][j]  + k14_10 * yDotK[9][j] +
          k14_11 * yDotK[10][j] + k14_12 * yDotK[11][j] + k14_13 * yDotK[12][j];
      yTmp[j] = currentState[j] + h * s;
    }
    equations.computeDerivatives(previousTime + c14 * h, yTmp, yDotKLast[0]);

    // k15
    for (int j = 0; j < currentState.length; ++j) {
     s = k15_01 * yDotK[0][j]  + k15_06 * yDotK[5][j]  + k15_07 * yDotK[6][j] +
         k15_08 * yDotK[7][j]  + k15_09 * yDotK[8][j]  + k15_10 * yDotK[9][j] +
         k15_11 * yDotK[10][j] + k15_12 * yDotK[11][j] + k15_13 * yDotK[12][j] +
         k15_14 * yDotKLast[0][j];
     yTmp[j] = currentState[j] + h * s;
    }
    equations.computeDerivatives(previousTime + c15 * h, yTmp, yDotKLast[1]);

    // k16
    for (int j = 0; j < currentState.length; ++j) {