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.ode.DerivativeException;
026    import org.apache.commons.math.linear.Array2DRowRealMatrix;
027    import org.apache.commons.math.util.FastMath;
028    
029    /**
030     * This class implements an interpolator for integrators using Nordsieck representation.
031     *
032     * <p>This interpolator computes dense output around the current point.
033     * The interpolation equation is based on Taylor series formulas.
034     *
035     * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator
036     * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator
037     * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 f??vr. 2011) $
038     * @since 2.0
039     */
040    
041    public 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         */
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    }