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