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
018package org.apache.commons.math4.legacy.ode.sampling;
019
020import java.io.IOException;
021import java.io.ObjectInput;
022import java.io.ObjectOutput;
023import java.util.Arrays;
024
025import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
026import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix;
027import org.apache.commons.math4.legacy.ode.EquationsMapper;
028import org.apache.commons.math4.core.jdkmath.JdkMath;
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.math4.legacy.ode.nonstiff.AdamsBashforthIntegrator
037 * @see org.apache.commons.math4.legacy.ode.nonstiff.AdamsMoultonIntegrator
038 * @since 2.0
039 */
040
041public 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     * @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    }
118
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;
136
137        // make sure the state and derivatives will depend on the new arrays
138        setInterpolatedTime(getInterpolatedTime());
139    }
140
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) {
147
148        final double ratio = stepSize / scalingH;
149        for (int i = 0; i < scaled.length; ++i) {
150            scaled[i] *= ratio;
151        }
152
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        }
162
163        scalingH = stepSize;
164    }
165
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    }
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 = 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        }
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    /** {@inheritDoc} */
218    @Override
219    public void writeExternal(final ObjectOutput out)
220        throws IOException {
221
222        // save the state of the base class
223        writeBaseExternal(out);
224
225        // save the local attributes
226        out.writeDouble(scalingH);
227        out.writeDouble(referenceTime);
228
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        }
238
239        if (nordsieck == null) {
240            out.writeBoolean(false);
241        } else {
242            out.writeBoolean(true);
243            out.writeObject(nordsieck);
244        }
245
246        // we don't save state variation, it will be recomputed
247    }
248
249    /** {@inheritDoc} */
250    @Override
251    public void readExternal(final ObjectInput in)
252        throws IOException, ClassNotFoundException {
253
254        // read the base class
255        final double t = readBaseExternal(in);
256
257        // read the local attributes
258        scalingH      = in.readDouble();
259        referenceTime = in.readDouble();
260
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        }
271
272        final boolean hasNordsieck = in.readBoolean();
273        if (hasNordsieck) {
274            nordsieck = (Array2DRowRealMatrix) in.readObject();
275        } else {
276            nordsieck = null;
277        }
278
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}