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.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.exception.MaxCountExceededException;
026import org.apache.commons.math4.linear.Array2DRowRealMatrix;
027import org.apache.commons.math4.ode.EquationsMapper;
028import org.apache.commons.math4.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.math4.ode.nonstiff.AdamsBashforthIntegrator
037 * @see org.apache.commons.math4.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
142    /** Rescale the instance.
143     * <p>Since the scaled and Nordsieck arrays are shared with the caller,
144     * this method has the side effect of rescaling this arrays in the caller too.</p>
145     * @param stepSize new step size to use in the scaled and Nordsieck arrays
146     */
147    public void rescale(final double stepSize) {
148
149        final double ratio = stepSize / scalingH;
150        for (int i = 0; i < scaled.length; ++i) {
151            scaled[i] *= ratio;
152        }
153
154        final double[][] nData = nordsieck.getDataRef();
155        double power = ratio;
156        for (int i = 0; i < nData.length; ++i) {
157            power *= ratio;
158            final double[] nDataI = nData[i];
159            for (int j = 0; j < nDataI.length; ++j) {
160                nDataI[j] *= power;
161            }
162        }
163
164        scalingH = stepSize;
165
166    }
167
168    /**
169     * Get the state vector variation from current to interpolated state.
170     * <p>This method is aimed at computing y(t<sub>interpolation</sub>)
171     * -y(t<sub>current</sub>) accurately by avoiding the cancellation errors
172     * that would occur if the subtraction were performed explicitly.</p>
173     * <p>The returned vector is a reference to a reused array, so
174     * it should not be modified and it should be copied if it needs
175     * to be preserved across several calls.</p>
176     * @return state vector at time {@link #getInterpolatedTime}
177     * @see #getInterpolatedDerivatives()
178     * @exception MaxCountExceededException if the number of functions evaluations is exceeded
179     */
180    public double[] getInterpolatedStateVariation() throws MaxCountExceededException {
181        // compute and ignore interpolated state
182        // to make sure state variation is computed as a side effect
183        getInterpolatedState();
184        return stateVariation;
185    }
186
187    /** {@inheritDoc} */
188    @Override
189    protected void computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH) {
190
191        final double x = interpolatedTime - referenceTime;
192        final double normalizedAbscissa = x / scalingH;
193
194        Arrays.fill(stateVariation, 0.0);
195        Arrays.fill(interpolatedDerivatives, 0.0);
196
197        // apply Taylor formula from high order to low order,
198        // for the sake of numerical accuracy
199        final double[][] nData = nordsieck.getDataRef();
200        for (int i = nData.length - 1; i >= 0; --i) {
201            final int order = i + 2;
202            final double[] nDataI = nData[i];
203            final double power = FastMath.pow(normalizedAbscissa, order);
204            for (int j = 0; j < nDataI.length; ++j) {
205                final double d = nDataI[j] * power;
206                stateVariation[j]          += d;
207                interpolatedDerivatives[j] += order * d;
208            }
209        }
210
211        for (int j = 0; j < currentState.length; ++j) {
212            stateVariation[j] += scaled[j] * normalizedAbscissa;
213            interpolatedState[j] = currentState[j] + stateVariation[j];
214            interpolatedDerivatives[j] =
215                (interpolatedDerivatives[j] + scaled[j] * normalizedAbscissa) / x;
216        }
217
218    }
219
220    /** {@inheritDoc} */
221    @Override
222    public void writeExternal(final ObjectOutput out)
223        throws IOException {
224
225        // save the state of the base class
226        writeBaseExternal(out);
227
228        // save the local attributes
229        out.writeDouble(scalingH);
230        out.writeDouble(referenceTime);
231
232        final int n = (currentState == null) ? -1 : currentState.length;
233        if (scaled == null) {
234            out.writeBoolean(false);
235        } else {
236            out.writeBoolean(true);
237            for (int j = 0; j < n; ++j) {
238                out.writeDouble(scaled[j]);
239            }
240        }
241
242        if (nordsieck == null) {
243            out.writeBoolean(false);
244        } else {
245            out.writeBoolean(true);
246            out.writeObject(nordsieck);
247        }
248
249        // we don't save state variation, it will be recomputed
250
251    }
252
253    /** {@inheritDoc} */
254    @Override
255    public void readExternal(final ObjectInput in)
256        throws IOException, ClassNotFoundException {
257
258        // read the base class
259        final double t = readBaseExternal(in);
260
261        // read the local attributes
262        scalingH      = in.readDouble();
263        referenceTime = in.readDouble();
264
265        final int n = (currentState == null) ? -1 : currentState.length;
266        final boolean hasScaled = in.readBoolean();
267        if (hasScaled) {
268            scaled = new double[n];
269            for (int j = 0; j < n; ++j) {
270                scaled[j] = in.readDouble();
271            }
272        } else {
273            scaled = null;
274        }
275
276        final boolean hasNordsieck = in.readBoolean();
277        if (hasNordsieck) {
278            nordsieck = (Array2DRowRealMatrix) in.readObject();
279        } else {
280            nordsieck = null;
281        }
282
283        if (hasScaled && hasNordsieck) {
284            // we can now set the interpolated time and state
285            stateVariation = new double[n];
286            setInterpolatedTime(t);
287        } else {
288            stateVariation = null;
289        }
290
291    }
292
293}