RungeKuttaStepInterpolator.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  17. package org.apache.commons.math4.legacy.ode.nonstiff;

  18. import java.io.IOException;
  19. import java.io.ObjectInput;
  20. import java.io.ObjectOutput;

  21. import org.apache.commons.math4.legacy.ode.AbstractIntegrator;
  22. import org.apache.commons.math4.legacy.ode.EquationsMapper;
  23. import org.apache.commons.math4.legacy.ode.sampling.AbstractStepInterpolator;

  24. /** This class represents an interpolator over the last step during an
  25.  * ODE integration for Runge-Kutta and embedded Runge-Kutta integrators.
  26.  *
  27.  * @see RungeKuttaIntegrator
  28.  * @see EmbeddedRungeKuttaIntegrator
  29.  *
  30.  * @since 1.2
  31.  */

  32. abstract class RungeKuttaStepInterpolator
  33.   extends AbstractStepInterpolator {

  34.     /** Previous state. */
  35.     protected double[] previousState;

  36.     /** Slopes at the intermediate points. */
  37.     protected double[][] yDotK;

  38.     /** Reference to the integrator. */
  39.     protected AbstractIntegrator integrator;

  40.   /** Simple constructor.
  41.    * This constructor builds an instance that is not usable yet, the
  42.    * {@link #reinitialize} method should be called before using the
  43.    * instance in order to initialize the internal arrays. This
  44.    * constructor is used only in order to delay the initialization in
  45.    * some cases. The {@link RungeKuttaIntegrator} and {@link
  46.    * EmbeddedRungeKuttaIntegrator} classes use the prototyping design
  47.    * pattern to create the step interpolators by cloning an
  48.    * uninitialized model and latter initializing the copy.
  49.    */
  50.   protected RungeKuttaStepInterpolator() {
  51.     previousState = null;
  52.     yDotK         = null;
  53.     integrator    = null;
  54.   }

  55.   /** Copy constructor.

  56.   * <p>The copied interpolator should have been finalized before the
  57.   * copy, otherwise the copy will not be able to perform correctly any
  58.   * interpolation and will throw a {@link NullPointerException}
  59.   * later. Since we don't want this constructor to throw the
  60.   * exceptions finalization may involve and since we don't want this
  61.   * method to modify the state of the copied interpolator,
  62.   * finalization is <strong>not</strong> done automatically, it
  63.   * remains under user control.</p>

  64.   * <p>The copy is a deep copy: its arrays are separated from the
  65.   * original arrays of the instance.</p>

  66.   * @param interpolator interpolator to copy from.

  67.   */
  68.   RungeKuttaStepInterpolator(final RungeKuttaStepInterpolator interpolator) {

  69.     super(interpolator);

  70.     if (interpolator.currentState != null) {

  71.       previousState = interpolator.previousState.clone();

  72.       yDotK = new double[interpolator.yDotK.length][];
  73.       for (int k = 0; k < interpolator.yDotK.length; ++k) {
  74.         yDotK[k] = interpolator.yDotK[k].clone();
  75.       }
  76.     } else {
  77.       previousState = null;
  78.       yDotK = null;
  79.     }

  80.     // we cannot keep any reference to the equations in the copy
  81.     // the interpolator should have been finalized before
  82.     integrator = null;
  83.   }

  84.   /** Reinitialize the instance
  85.    * <p>Some Runge-Kutta integrators need fewer functions evaluations
  86.    * than their counterpart step interpolators. So the interpolator
  87.    * should perform the last evaluations they need by themselves. The
  88.    * {@link RungeKuttaIntegrator RungeKuttaIntegrator} and {@link
  89.    * EmbeddedRungeKuttaIntegrator EmbeddedRungeKuttaIntegrator}
  90.    * abstract classes call this method in order to let the step
  91.    * interpolator perform the evaluations it needs. These evaluations
  92.    * will be performed during the call to <code>doFinalize</code> if
  93.    * any, i.e. only if the step handler either calls the {@link
  94.    * AbstractStepInterpolator#finalizeStep finalizeStep} method or the
  95.    * {@link AbstractStepInterpolator#getInterpolatedState
  96.    * getInterpolatedState} method (for an interpolator which needs a
  97.    * finalization) or if it clones the step interpolator.</p>
  98.    * @param rkIntegrator integrator being used
  99.    * @param y reference to the integrator array holding the state at
  100.    * the end of the step
  101.    * @param yDotArray reference to the integrator array holding all the
  102.    * intermediate slopes
  103.    * @param forward integration direction indicator
  104.    * @param primaryMapper equations mapper for the primary equations set
  105.    * @param secondaryMappers equations mappers for the secondary equations sets
  106.    */
  107.   public void reinitialize(final AbstractIntegrator rkIntegrator,
  108.                            final double[] y, final double[][] yDotArray, final boolean forward,
  109.                            final EquationsMapper primaryMapper,
  110.                            final EquationsMapper[] secondaryMappers) {
  111.     reinitialize(y, forward, primaryMapper, secondaryMappers);
  112.     this.previousState = null;
  113.     this.yDotK = yDotArray;
  114.     this.integrator = rkIntegrator;
  115.   }

  116.   /** {@inheritDoc} */
  117.   @Override
  118.   public void shift() {
  119.     previousState = currentState.clone();
  120.     super.shift();
  121.   }

  122.   /** {@inheritDoc} */
  123.   @Override
  124.   public void writeExternal(final ObjectOutput out)
  125.     throws IOException {

  126.     // save the state of the base class
  127.     writeBaseExternal(out);

  128.     // save the local attributes
  129.     final int n = (currentState == null) ? -1 : currentState.length;
  130.     for (int i = 0; i < n; ++i) {
  131.       out.writeDouble(previousState[i]);
  132.     }

  133.     final int kMax = (yDotK == null) ? -1 : yDotK.length;
  134.     out.writeInt(kMax);
  135.     for (int k = 0; k < kMax; ++k) {
  136.       for (int i = 0; i < n; ++i) {
  137.         out.writeDouble(yDotK[k][i]);
  138.       }
  139.     }

  140.     // we do not save any reference to the equations
  141.   }

  142.   /** {@inheritDoc} */
  143.   @Override
  144.   public void readExternal(final ObjectInput in)
  145.     throws IOException, ClassNotFoundException {

  146.     // read the base class
  147.     final double t = readBaseExternal(in);

  148.     // read the local attributes
  149.     final int n = (currentState == null) ? -1 : currentState.length;
  150.     if (n < 0) {
  151.       previousState = null;
  152.     } else {
  153.       previousState = new double[n];
  154.       for (int i = 0; i < n; ++i) {
  155.         previousState[i] = in.readDouble();
  156.       }
  157.     }

  158.     final int kMax = in.readInt();
  159.     yDotK = (kMax < 0) ? null : new double[kMax][];
  160.     for (int k = 0; k < kMax; ++k) {
  161.       yDotK[k] = (n < 0) ? null : new double[n];
  162.       for (int i = 0; i < n; ++i) {
  163.         yDotK[k][i] = in.readDouble();
  164.       }
  165.     }

  166.     integrator = null;

  167.     if (currentState != null) {
  168.         // we can now set the interpolated time and state
  169.         setInterpolatedTime(t);
  170.     } else {
  171.         interpolatedTime = t;
  172.     }
  173.   }
  174. }