    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017
    018    package org.apache.commons.math3.ode.sampling;
    019
    020    import java.io.IOException;
    021    import java.io.ObjectInput;
    022    import java.io.ObjectOutput;
    023    import java.util.Arrays;
    024
    025    import org.apache.commons.math3.exception.MaxCountExceededException;
    026    import org.apache.commons.math3.linear.Array2DRowRealMatrix;
    027    import org.apache.commons.math3.ode.EquationsMapper;
    028    import org.apache.commons.math3.util.FastMath;
    029
    030    /**
    031     * This class implements an interpolator for integrators using Nordsieck representation.
    032     *
    033     * <p>This interpolator computes dense output around the current point.
    034     * The interpolation equation is based on Taylor series formulas.
    035     *
    036     * @see org.apache.commons.math3.ode.nonstiff.AdamsBashforthIntegrator
    037     * @see org.apache.commons.math3.ode.nonstiff.AdamsMoultonIntegrator
    038     * @version $Id: NordsieckStepInterpolator.java 1416643 2012-12-03 19:37:14Z tn $
    039     * @since 2.0
    040     */
    041
    042    public class NordsieckStepInterpolator extends AbstractStepInterpolator {
    043
    044        /** Serializable version identifier */
    045        private static final long serialVersionUID = -7179861704951334960L;
    046
    047        /** State variation. */
    048        protected double[] stateVariation;
    049
    050        /** Step size used in the first scaled derivative and Nordsieck vector. */
    051        private double scalingH;
    052
    053        /** Reference time for all arrays.
    054         * <p>Sometimes, the reference time is the same as previousTime,
    055         * sometimes it is the same as currentTime, so we use a separate
    056         * field to avoid any confusion.
    057         * </p>
    058         */
    059        private double referenceTime;
    060
    061        /** First scaled derivative. */
    062        private double[] scaled;
    063
    064        /** Nordsieck vector. */
    065        private Array2DRowRealMatrix nordsieck;
    066
    067        /** Simple constructor.
    068         * This constructor builds an instance that is not usable yet, the
    069         * {@link AbstractStepInterpolator#reinitialize} method should be called
    070         * before using the instance in order to initialize the internal arrays. This
    071         * constructor is used only in order to delay the initialization in
    072         * some cases.
    073         */
    074        public NordsieckStepInterpolator() {
    075        }
    076
    077        /** Copy constructor.
    078         * @param interpolator interpolator to copy from. The copy is a deep
    079         * copy: its arrays are separated from the original arrays of the
    080         * instance
    081         */
    082        public NordsieckStepInterpolator(final NordsieckStepInterpolator interpolator) {
    083            super(interpolator);
    084            scalingH      = interpolator.scalingH;
    085            referenceTime = interpolator.referenceTime;
    086            if (interpolator.scaled != null) {
    087                scaled = interpolator.scaled.clone();
    088            }
    089            if (interpolator.nordsieck != null) {
    090                nordsieck = new Array2DRowRealMatrix(interpolator.nordsieck.getDataRef(), true);
    091            }
    092            if (interpolator.stateVariation != null) {
    093                stateVariation = interpolator.stateVariation.clone();
    094            }
    095        }
    096
    097        /** {@inheritDoc} */
    098        @Override
    099        protected StepInterpolator doCopy() {
    100            return new NordsieckStepInterpolator(this);
    101        }
    102
    103        /** Reinitialize the instance.
    104         * <p>Beware that all arrays <em>must</em> be references to integrator
    105         * arrays, in order to ensure proper update without copy.</p>
    106         * @param y reference to the integrator array holding the state at
    107         * the end of the step
    108         * @param forward integration direction indicator
    109         * @param primaryMapper equations mapper for the primary equations set
    110         * @param secondaryMappers equations mappers for the secondary equations sets
    111         */
    112        @Override
    113        public void reinitialize(final double[] y, final boolean forward,
    114                                 final EquationsMapper primaryMapper,
    115                                 final EquationsMapper[] secondaryMappers) {
    116            super.reinitialize(y, forward, primaryMapper, secondaryMappers);
    117            stateVariation = new double[y.length];
    118        }
    119
    120        /** Reinitialize the instance.
    121         * <p>Beware that all arrays <em>must</em> be references to integrator
    122         * arrays, in order to ensure proper update without copy.</p>
    123         * @param time time at which all arrays are defined
    124         * @param stepSize step size used in the scaled and nordsieck arrays
    125         * @param scaledDerivative reference to the integrator array holding the first
    126         * scaled derivative
    127         * @param nordsieckVector reference to the integrator matrix holding the
    128         * nordsieck vector
    129         */
    130        public void reinitialize(final double time, final double stepSize,
    131                                 final double[] scaledDerivative,
    132                                 final Array2DRowRealMatrix nordsieckVector) {
    133            this.referenceTime = time;
    134            this.scalingH      = stepSize;
    135            this.scaled        = scaledDerivative;
    136            this.nordsieck     = nordsieckVector;
    137
    138            // make sure the state and derivatives will depend on the new arrays
    139            setInterpolatedTime(getInterpolatedTime());
    140
    141        }
    142
    143        /** Rescale the instance.
    144         * <p>Since the scaled and Nordiseck arrays are shared with the caller,
    145         * this method has the side effect of rescaling this arrays in the caller too.</p>
    146         * @param stepSize new step size to use in the scaled and nordsieck arrays
    147         */
    148        public void rescale(final double stepSize) {
    149
    150            final double ratio = stepSize / scalingH;
    151            for (int i = 0; i < scaled.length; ++i) {
    152                scaled[i] *= ratio;
    153            }
    154
    155            final double[][] nData = nordsieck.getDataRef();
    156            double power = ratio;
    157            for (int i = 0; i < nData.length; ++i) {
    158                power *= ratio;
    159                final double[] nDataI = nData[i];
    160                for (int j = 0; j < nDataI.length; ++j) {
    161                    nDataI[j] *= power;
    162                }
    163            }
    164
    165            scalingH = stepSize;
    166
    167        }
    168
    169        /**
    170         * Get the state vector variation from current to interpolated state.
    171         * <p>This method is aimed at computing y(t<sub>interpolation</sub>)
    172         * -y(t<sub>current</sub>) accurately by avoiding the cancellation errors
    173         * that would occur if the subtraction were performed explicitly.</p>
    174         * <p>The returned vector is a reference to a reused array, so
    175         * it should not be modified and it should be copied if it needs
    176         * to be preserved across several calls.</p>
    177         * @return state vector at time {@link #getInterpolatedTime}
    178         * @see #getInterpolatedDerivatives()
    179         * @exception MaxCountExceededException if the number of functions evaluations is exceeded
    180         */
    181        public double[] getInterpolatedStateVariation() throws MaxCountExceededException {
    182            // compute and ignore interpolated state
    183            // to make sure state variation is computed as a side effect
    184            getInterpolatedState();
    185            return stateVariation;
    186        }
    187
    188        /** {@inheritDoc} */
    189        @Override
    190        protected void computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH) {
    191
    192            final double x = interpolatedTime - referenceTime;
    193            final double normalizedAbscissa = x / scalingH;
    194
    195            Arrays.fill(stateVariation, 0.0);
    196            Arrays.fill(interpolatedDerivatives, 0.0);
    197
    198            // apply Taylor formula from high order to low order,
    199            // for the sake of numerical accuracy
    200            final double[][] nData = nordsieck.getDataRef();
    201            for (int i = nData.length - 1; i >= 0; --i) {
    202                final int order = i + 2;
    203                final double[] nDataI = nData[i];
    204                final double power = FastMath.pow(normalizedAbscissa, order);
    205                for (int j = 0; j < nDataI.length; ++j) {
    206                    final double d = nDataI[j] * power;
    207                    stateVariation[j]          += d;
    208                    interpolatedDerivatives[j] += order * d;
    209                }
    210            }
    211
    212            for (int j = 0; j < currentState.length; ++j) {
    213                stateVariation[j] += scaled[j] * normalizedAbscissa;
    214                interpolatedState[j] = currentState[j] + stateVariation[j];
    215                interpolatedDerivatives[j] =
    216                    (interpolatedDerivatives[j] + scaled[j] * normalizedAbscissa) / x;
    217            }
    218
    219        }
    220
    221        /** {@inheritDoc} */
    222        @Override
    223        public void writeExternal(final ObjectOutput out)
    224            throws IOException {
    225
    226            // save the state of the base class
    227            writeBaseExternal(out);
    228
    229            // save the local attributes
    230            out.writeDouble(scalingH);
    231            out.writeDouble(referenceTime);
    232
    233            final int n = (currentState == null) ? -1 : currentState.length;
    234            if (scaled == null) {
    235                out.writeBoolean(false);
    236            } else {
    237                out.writeBoolean(true);
    238                for (int j = 0; j < n; ++j) {
    239                    out.writeDouble(scaled[j]);
    240                }
    241            }
    242
    243            if (nordsieck == null) {
    244                out.writeBoolean(false);
    245            } else {
    246                out.writeBoolean(true);
    247                out.writeObject(nordsieck);
    248            }
    249
    250            // we don't save state variation, it will be recomputed
    251
    252        }
    253
    254        /** {@inheritDoc} */
    255        @Override
    256        public void readExternal(final ObjectInput in)
    257            throws IOException, ClassNotFoundException {
    258
    259            // read the base class
    260            final double t = readBaseExternal(in);
    261
    262            // read the local attributes
    263            scalingH      = in.readDouble();
    264            referenceTime = in.readDouble();
    265
    266            final int n = (currentState == null) ? -1 : currentState.length;
    267            final boolean hasScaled = in.readBoolean();
    268            if (hasScaled) {
    269                scaled = new double[n];
    270                for (int j = 0; j < n; ++j) {
    271                    scaled[j] = in.readDouble();
    272                }
    273            } else {
    274                scaled = null;
    275            }
    276
    277            final boolean hasNordsieck = in.readBoolean();
    278            if (hasNordsieck) {
    279                nordsieck = (Array2DRowRealMatrix) in.readObject();
    280            } else {
    281                nordsieck = null;
    282            }
    283
    284            if (hasScaled && hasNordsieck) {
    285                // we can now set the interpolated time and state
    286                stateVariation = new double[n];
    287                setInterpolatedTime(t);
    288            } else {
    289                stateVariation = null;
    290            }
    291
    292        }
    293
    294    }