1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements.  See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License.  You may obtain a copy of the License at
8 *
9 *      http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18package org.apache.commons.math.ode.sampling;
19
20import java.io.IOException;
21import java.io.ObjectInput;
22import java.io.ObjectOutput;
23import java.util.Arrays;
24
25import org.apache.commons.math.ode.DerivativeException;
26import org.apache.commons.math.linear.Array2DRowRealMatrix;
27import org.apache.commons.math.util.FastMath;
28
29/**
30 * This class implements an interpolator for integrators using Nordsieck representation.
31 *
32 * <p>This interpolator computes dense output around the current point.
33 * The interpolation equation is based on Taylor series formulas.
34 *
35 * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator
36 * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator
37 * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $
38 * @since 2.0
39 */
40
41public class NordsieckStepInterpolator extends AbstractStepInterpolator {
42
43    /** Serializable version identifier */
44    private static final long serialVersionUID = -7179861704951334960L;
45
46    /** State variation. */
47    protected double[] stateVariation;
48
49    /** Step size used in the first scaled derivative and Nordsieck vector. */
50    private double scalingH;
51
52    /** Reference time for all arrays.
53     * <p>Sometimes, the reference time is the same as previousTime,
54     * sometimes it is the same as currentTime, so we use a separate
55     * field to avoid any confusion.
56     * </p>
57     */
58    private double referenceTime;
59
60    /** First scaled derivative. */
61    private double[] scaled;
62
63    /** Nordsieck vector. */
64    private Array2DRowRealMatrix nordsieck;
65
66    /** Simple constructor.
67     * This constructor builds an instance that is not usable yet, the
68     * {@link AbstractStepInterpolator#reinitialize} method should be called
69     * before using the instance in order to initialize the internal arrays. This
70     * constructor is used only in order to delay the initialization in
71     * some cases.
72     */
73    public NordsieckStepInterpolator() {
74    }
75
76    /** Copy constructor.
77     * @param interpolator interpolator to copy from. The copy is a deep
78     * copy: its arrays are separated from the original arrays of the
79     * instance
80     */
81    public NordsieckStepInterpolator(final NordsieckStepInterpolator interpolator) {
82        super(interpolator);
83        scalingH      = interpolator.scalingH;
84        referenceTime = interpolator.referenceTime;
85        if (interpolator.scaled != null) {
86            scaled = interpolator.scaled.clone();
87        }
88        if (interpolator.nordsieck != null) {
89            nordsieck = new Array2DRowRealMatrix(interpolator.nordsieck.getDataRef(), true);
90        }
91        if (interpolator.stateVariation != null) {
92            stateVariation = interpolator.stateVariation.clone();
93        }
94    }
95
96    /** {@inheritDoc} */
97    @Override
98    protected StepInterpolator doCopy() {
99        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}
292