NordsieckStepInterpolator.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.  */

  18. package org.apache.commons.math4.legacy.ode.sampling;

  20. import java.io.IOException;
  21. import java.io.ObjectInput;
  22. import java.io.ObjectOutput;
  23. import java.util.Arrays;

  25. import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
  26. import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix;
  27. import org.apache.commons.math4.legacy.ode.EquationsMapper;
  28. import org.apache.commons.math4.core.jdkmath.JdkMath;

  30. /**
  31.  * This class implements an interpolator for integrators using Nordsieck representation.
  32.  *
  33.  * <p>This interpolator computes dense output around the current point.
  34.  * The interpolation equation is based on Taylor series formulas.
  35.  *
  36.  * @see org.apache.commons.math4.legacy.ode.nonstiff.AdamsBashforthIntegrator
  37.  * @see org.apache.commons.math4.legacy.ode.nonstiff.AdamsMoultonIntegrator
  38.  * @since 2.0
  39.  */

  41. public class NordsieckStepInterpolator extends AbstractStepInterpolator {

  43.     /** Serializable version identifier. */
  44.     private static final long serialVersionUID = -7179861704951334960L;

  46.     /** State variation. */
  47.     protected double[] stateVariation;

  49.     /** Step size used in the first scaled derivative and Nordsieck vector. */
  50.     private double scalingH;

  52.     /** Reference time for all arrays.
  53.      * <p>Sometimes, the reference time is the same as previousTime,
  54.      * sometimes it is the same as currentTime, so we use a separate
  55.      * field to avoid any confusion.
  56.      * </p>
  57.      */
  58.     private double referenceTime;

  60.     /** First scaled derivative. */
  61.     private double[] scaled;

  63.     /** Nordsieck vector. */
  64.     private Array2DRowRealMatrix nordsieck;

  66.     /** Simple constructor.
  67.      * This constructor builds an instance that is not usable yet, the
  68.      * {@link AbstractStepInterpolator#reinitialize} method should be called
  69.      * before using the instance in order to initialize the internal arrays. This
  70.      * constructor is used only in order to delay the initialization in
  71.      * some cases.
  72.      */
  73.     public NordsieckStepInterpolator() {
  74.     }

  76.     /** Copy constructor.
  77.      * @param interpolator interpolator to copy from. The copy is a deep
  78.      * copy: its arrays are separated from the original arrays of the
  79.      * instance
  80.      */
  81.     public NordsieckStepInterpolator(final NordsieckStepInterpolator interpolator) {
  82.         super(interpolator);
  83.         scalingH      = interpolator.scalingH;
  84.         referenceTime = interpolator.referenceTime;
  85.         if (interpolator.scaled != null) {
  86.             scaled = interpolator.scaled.clone();
  87.         }
  88.         if (interpolator.nordsieck != null) {
  89.             nordsieck = new Array2DRowRealMatrix(interpolator.nordsieck.getDataRef(), true);
  90.         }
  91.         if (interpolator.stateVariation != null) {
  92.             stateVariation = interpolator.stateVariation.clone();
  93.         }
  94.     }

  96.     /** {@inheritDoc} */
  97.     @Override
  98.     protected StepInterpolator doCopy() {
  99.         return new NordsieckStepInterpolator(this);
  100.     }

  102.     /** Reinitialize the instance.
  103.      * <p>Beware that all arrays <em>must</em> be references to integrator
  104.      * arrays, in order to ensure proper update without copy.</p>
  105.      * @param y reference to the integrator array holding the state at
  106.      * the end of the step
  107.      * @param forward integration direction indicator
  108.      * @param primaryMapper equations mapper for the primary equations set
  109.      * @param secondaryMappers equations mappers for the secondary equations sets
  110.      */
  111.     @Override
  112.     public void reinitialize(final double[] y, final boolean forward,
  113.                              final EquationsMapper primaryMapper,
  114.                              final EquationsMapper[] secondaryMappers) {
  115.         super.reinitialize(y, forward, primaryMapper, secondaryMappers);
  116.         stateVariation = new double[y.length];
  117.     }

  119.     /** Reinitialize the instance.
  120.      * <p>Beware that all arrays <em>must</em> be references to integrator
  121.      * arrays, in order to ensure proper update without copy.</p>
  122.      * @param time time at which all arrays are defined
  123.      * @param stepSize step size used in the scaled and Nordsieck arrays
  124.      * @param scaledDerivative reference to the integrator array holding the first
  125.      * scaled derivative
  126.      * @param nordsieckVector reference to the integrator matrix holding the
  127.      * Nordsieck vector
  128.      */
  129.     public void reinitialize(final double time, final double stepSize,
  130.                              final double[] scaledDerivative,
  131.                              final Array2DRowRealMatrix nordsieckVector) {
  132.         this.referenceTime = time;
  133.         this.scalingH      = stepSize;
  134.         this.scaled        = scaledDerivative;
  135.         this.nordsieck     = nordsieckVector;

  137.         // make sure the state and derivatives will depend on the new arrays
  138.         setInterpolatedTime(getInterpolatedTime());
  139.     }

  141.     /** Rescale the instance.
  142.      * <p>Since the scaled and Nordsieck arrays are shared with the caller,
  143.      * this method has the side effect of rescaling this arrays in the caller too.</p>
  144.      * @param stepSize new step size to use in the scaled and Nordsieck arrays
  145.      */
  146.     public void rescale(final double stepSize) {

  148.         final double ratio = stepSize / scalingH;
  149.         for (int i = 0; i < scaled.length; ++i) {
  150.             scaled[i] *= ratio;
  151.         }

  153.         final double[][] nData = nordsieck.getDataRef();
  154.         double power = ratio;
  155.         for (int i = 0; i < nData.length; ++i) {
  156.             power *= ratio;
  157.             final double[] nDataI = nData[i];
  158.             for (int j = 0; j < nDataI.length; ++j) {
  159.                 nDataI[j] *= power;
  160.             }
  161.         }

  163.         scalingH = stepSize;
  164.     }

  166.     /**
  167.      * Get the state vector variation from current to interpolated state.
  168.      * <p>This method is aimed at computing y(t<sub>interpolation</sub>)
  169.      * -y(t<sub>current</sub>) accurately by avoiding the cancellation errors
  170.      * that would occur if the subtraction were performed explicitly.</p>
  171.      * <p>The returned vector is a reference to a reused array, so
  172.      * it should not be modified and it should be copied if it needs
  173.      * to be preserved across several calls.</p>
  174.      * @return state vector at time {@link #getInterpolatedTime}
  175.      * @see #getInterpolatedDerivatives()
  176.      * @exception MaxCountExceededException if the number of functions evaluations is exceeded
  177.      */
  178.     public double[] getInterpolatedStateVariation() throws MaxCountExceededException {
  179.         // compute and ignore interpolated state
  180.         // to make sure state variation is computed as a side effect
  181.         getInterpolatedState();
  182.         return stateVariation;
  183.     }

  185.     /** {@inheritDoc} */
  186.     @Override
  187.     protected void computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH) {

  189.         final double x = interpolatedTime - referenceTime;
  190.         final double normalizedAbscissa = x / scalingH;

  192.         Arrays.fill(stateVariation, 0.0);
  193.         Arrays.fill(interpolatedDerivatives, 0.0);

  195.         // apply Taylor formula from high order to low order,
  196.         // for the sake of numerical accuracy
  197.         final double[][] nData = nordsieck.getDataRef();
  198.         for (int i = nData.length - 1; i >= 0; --i) {
  199.             final int order = i + 2;
  200.             final double[] nDataI = nData[i];
  201.             final double power = JdkMath.pow(normalizedAbscissa, order);
  202.             for (int j = 0; j < nDataI.length; ++j) {
  203.                 final double d = nDataI[j] * power;
  204.                 stateVariation[j]          += d;
  205.                 interpolatedDerivatives[j] += order * d;
  206.             }
  207.         }

  209.         for (int j = 0; j < currentState.length; ++j) {
  210.             stateVariation[j] += scaled[j] * normalizedAbscissa;
  211.             interpolatedState[j] = currentState[j] + stateVariation[j];
  212.             interpolatedDerivatives[j] =
  213.                 (interpolatedDerivatives[j] + scaled[j] * normalizedAbscissa) / x;
  214.         }
  215.     }

  217.     /** {@inheritDoc} */
  218.     @Override
  219.     public void writeExternal(final ObjectOutput out)
  220.         throws IOException {

  222.         // save the state of the base class
  223.         writeBaseExternal(out);

  225.         // save the local attributes
  226.         out.writeDouble(scalingH);
  227.         out.writeDouble(referenceTime);

  229.         final int n = (currentState == null) ? -1 : currentState.length;
  230.         if (scaled == null) {
  231.             out.writeBoolean(false);
  232.         } else {
  233.             out.writeBoolean(true);
  234.             for (int j = 0; j < n; ++j) {
  235.                 out.writeDouble(scaled[j]);
  236.             }
  237.         }

  239.         if (nordsieck == null) {
  240.             out.writeBoolean(false);
  241.         } else {
  242.             out.writeBoolean(true);
  243.             out.writeObject(nordsieck);
  244.         }

  246.         // we don't save state variation, it will be recomputed
  247.     }

  249.     /** {@inheritDoc} */
  250.     @Override
  251.     public void readExternal(final ObjectInput in)
  252.         throws IOException, ClassNotFoundException {

  254.         // read the base class
  255.         final double t = readBaseExternal(in);

  257.         // read the local attributes
  258.         scalingH      = in.readDouble();
  259.         referenceTime = in.readDouble();

  261.         final int n = (currentState == null) ? -1 : currentState.length;
  262.         final boolean hasScaled = in.readBoolean();
  263.         if (hasScaled) {
  264.             scaled = new double[n];
  265.             for (int j = 0; j < n; ++j) {
  266.                 scaled[j] = in.readDouble();
  267.             }
  268.         } else {
  269.             scaled = null;
  270.         }

  272.         final boolean hasNordsieck = in.readBoolean();
  273.         if (hasNordsieck) {
  274.             nordsieck = (Array2DRowRealMatrix) in.readObject();
  275.         } else {
  276.             nordsieck = null;
  277.         }

  279.         if (hasScaled && hasNordsieck) {
  280.             // we can now set the interpolated time and state
  281.             stateVariation = new double[n];
  282.             setInterpolatedTime(t);
  283.         } else {
  284.             stateVariation = null;
  285.         }
  286.     }
  287. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT