    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.math.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.math.ode.DerivativeException;
    026    import org.apache.commons.math.linear.Array2DRowRealMatrix;
    027    import org.apache.commons.math.util.FastMath;
    028
    029    /**
    030     * This class implements an interpolator for integrators using Nordsieck representation.
    031     *
    032     * <p>This interpolator computes dense output around the current point.
    033     * The interpolation equation is based on Taylor series formulas.
    034     *
    035     * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator
    036     * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator
    037     * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 f??vr. 2011) $
    038     * @since 2.0
    039     */
    040
    041    public class NordsieckStepInterpolator extends AbstractStepInterpolator {
    042
    043        /** Serializable version identifier */
    044        private static final long serialVersionUID = -7179861704951334960L;
    045
    046        /** State variation. */
    047        protected double[] stateVariation;
    048
    049        /** Step size used in the first scaled derivative and Nordsieck vector. */
    050        private double scalingH;
    051
    052        /** Reference time for all arrays.
    053         * <p>Sometimes, the reference time is the same as previousTime,
    054         * sometimes it is the same as currentTime, so we use a separate
    055         * field to avoid any confusion.
    056         * </p>
    057         */
    058        private double referenceTime;
    059
    060        /** First scaled derivative. */
    061        private double[] scaled;
    062
    063        /** Nordsieck vector. */
    064        private Array2DRowRealMatrix nordsieck;
    065
    066        /** Simple constructor.
    067         * This constructor builds an instance that is not usable yet, the
    068         * {@link AbstractStepInterpolator#reinitialize} method should be called
    069         * before using the instance in order to initialize the internal arrays. This
    070         * constructor is used only in order to delay the initialization in
    071         * some cases.
    072         */
    073        public NordsieckStepInterpolator() {
    074        }
    075
    076        /** Copy constructor.
    077         * @param interpolator interpolator to copy from. The copy is a deep
    078         * copy: its arrays are separated from the original arrays of the
    079         * instance
    080         */
    081        public NordsieckStepInterpolator(final NordsieckStepInterpolator interpolator) {
    082            super(interpolator);
    083            scalingH      = interpolator.scalingH;
    084            referenceTime = interpolator.referenceTime;
    085            if (interpolator.scaled != null) {
    086                scaled = interpolator.scaled.clone();
    087            }
    088            if (interpolator.nordsieck != null) {
    089                nordsieck = new Array2DRowRealMatrix(interpolator.nordsieck.getDataRef(), true);
    090            }
    091            if (interpolator.stateVariation != null) {
    092                stateVariation = interpolator.stateVariation.clone();
    093            }
    094        }
    095
    096        /** {@inheritDoc} */
    097        @Override
    098        protected StepInterpolator doCopy() {
    099            return new NordsieckStepInterpolator(this);
    100        }
    101
    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         */
    109        @Override
    110        public void reinitialize(final double[] y, final boolean forward) {
    111            super.reinitialize(y, forward);
    112            stateVariation = new double[y.length];
    113        }
    114
    115        /** Reinitialize the instance.
    116         * <p>Beware that all arrays <em>must</em> be references to integrator
    117         * arrays, in order to ensure proper update without copy.</p>
    118         * @param time time at which all arrays are defined
    119         * @param stepSize step size used in the scaled and nordsieck arrays
    120         * @param scaledDerivative reference to the integrator array holding the first
    121         * scaled derivative
    122         * @param nordsieckVector reference to the integrator matrix holding the
    123         * nordsieck vector
    124         */
    125        public void reinitialize(final double time, final double stepSize,
    126                                 final double[] scaledDerivative,
    127                                 final Array2DRowRealMatrix nordsieckVector) {
    128            this.referenceTime = time;
    129            this.scalingH      = stepSize;
    130            this.scaled        = scaledDerivative;
    131            this.nordsieck     = nordsieckVector;
    132
    133            // make sure the state and derivatives will depend on the new arrays
    134            setInterpolatedTime(getInterpolatedTime());
    135
    136        }
    137
    138        /** Rescale the instance.
    139         * <p>Since the scaled and Nordiseck arrays are shared with the caller,
    140         * this method has the side effect of rescaling this arrays in the caller too.</p>
    141         * @param stepSize new step size to use in the scaled and nordsieck arrays
    142         */
    143        public void rescale(final double stepSize) {
    144
    145            final double ratio = stepSize / scalingH;
    146            for (int i = 0; i < scaled.length; ++i) {
    147                scaled[i] *= ratio;
    148            }
    149
    150            final double[][] nData = nordsieck.getDataRef();
    151            double power = ratio;
    152            for (int i = 0; i < nData.length; ++i) {
    153                power *= ratio;
    154                final double[] nDataI = nData[i];
    155                for (int j = 0; j < nDataI.length; ++j) {
    156                    nDataI[j] *= power;
    157                }
    158            }
    159
    160            scalingH = stepSize;
    161
    162        }
    163
    164        /**
    165         * Get the state vector variation from current to interpolated state.
    166         * <p>This method is aimed at computing y(t<sub>interpolation</sub>)
    167         * -y(t<sub>current</sub>) accurately by avoiding the cancellation errors
    168         * that would occur if the subtraction were performed explicitly.</p>
    169         * <p>The returned vector is a reference to a reused array, so
    170         * it should not be modified and it should be copied if it needs
    171         * to be preserved across several calls.</p>
    172         * @return state vector at time {@link #getInterpolatedTime}
    173         * @see #getInterpolatedDerivatives()
    174         * @throws DerivativeException if this call induces an automatic
    175         * step finalization that throws one
    176         */
    177        public double[] getInterpolatedStateVariation()
    178            throws DerivativeException {
    179            // compute and ignore interpolated state
    180            // to make sure state variation is computed as a side effect
    181            getInterpolatedState();
    182            return stateVariation;
    183        }
    184
    185        /** {@inheritDoc} */
    186        @Override
    187        protected void computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH) {
    188
    189            final double x = interpolatedTime - referenceTime;
    190            final double normalizedAbscissa = x / scalingH;
    191
    192            Arrays.fill(stateVariation, 0.0);
    193            Arrays.fill(interpolatedDerivatives, 0.0);
    194
    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 = FastMath.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            }
    208
    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
    216        }
    217
    218        /** {@inheritDoc} */
    219        @Override
    220        public void writeExternal(final ObjectOutput out)
    221            throws IOException {
    222
    223            // save the state of the base class
    224            writeBaseExternal(out);
    225
    226            // save the local attributes
    227            out.writeDouble(scalingH);
    228            out.writeDouble(referenceTime);
    229
    230            final int n = (currentState == null) ? -1 : currentState.length;
    231            if (scaled == null) {
    232                out.writeBoolean(false);
    233            } else {
    234                out.writeBoolean(true);
    235                for (int j = 0; j < n; ++j) {
    236                    out.writeDouble(scaled[j]);
    237                }
    238            }
    239
    240            if (nordsieck == null) {
    241                out.writeBoolean(false);
    242            } else {
    243                out.writeBoolean(true);
    244                out.writeObject(nordsieck);
    245            }
    246
    247            // we don't save state variation, it will be recomputed
    248
    249        }
    250
    251        /** {@inheritDoc} */
    252        @Override
    253        public void readExternal(final ObjectInput in)
    254            throws IOException, ClassNotFoundException {
    255
    256            // read the base class
    257            final double t = readBaseExternal(in);
    258
    259            // read the local attributes
    260            scalingH      = in.readDouble();
    261            referenceTime = in.readDouble();
    262
    263            final int n = (currentState == null) ? -1 : currentState.length;
    264            final boolean hasScaled = in.readBoolean();
    265            if (hasScaled) {
    266                scaled = new double[n];
    267                for (int j = 0; j < n; ++j) {
    268                    scaled[j] = in.readDouble();
    269                }
    270            } else {
    271                scaled = null;
    272            }
    273
    274            final boolean hasNordsieck = in.readBoolean();
    275            if (hasNordsieck) {
    276                nordsieck = (Array2DRowRealMatrix) in.readObject();
    277            } else {
    278                nordsieck = null;
    279            }
    280
    281            if (hasScaled && hasNordsieck) {
    282                // we can now set the interpolated time and state
    283                stateVariation = new double[n];
    284                setInterpolatedTime(t);
    285            } else {
    286                stateVariation = null;
    287            }
    288
    289        }
    290
    291    }