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.math3.ode.sampling;
019
020import java.io.IOException;
021import java.io.ObjectInput;
022import java.io.ObjectOutput;
023import java.util.Arrays;
024
025import org.apache.commons.math3.exception.MaxCountExceededException;
026import org.apache.commons.math3.linear.Array2DRowRealMatrix;
027import org.apache.commons.math3.ode.EquationsMapper;
028import 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 1503712 2013-07-16 13:35:19Z luc $
039 * @since 2.0
040 */
041
042public 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 Nordsieck 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}