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: 791723 $ $Date: 2009-07-07 03:07:23 -0400 (Tue, 07 Jul 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        /** Step size used in the first scaled derivative and Nordsieck vector. */
046        private double scalingH;
047    
048        /** Reference time for all arrays.
049         * <p>Sometimes, the reference time is the same as previousTime,
050         * sometimes it is the same as currentTime, so we use a separate
051         * field to avoid any confusion.
052         * </p>
053         */
054        private double referenceTime;
055    
056        /** First scaled derivative. */
057        private double[] scaled;
058    
059        /** Nordsieck vector. */
060        private Array2DRowRealMatrix nordsieck;
061    
062        /** State variation. */
063        protected double[] stateVariation;
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 referenceTime time at which all arrays are defined
118         * @param scalingH step size used in the scaled and nordsieck arrays
119         * @param scaled reference to the integrator array holding the first
120         * scaled derivative
121         * @param nordsieck reference to the integrator matrix holding the
122         * nordsieck vector
123         */
124        public void reinitialize(final double referenceTime, final double scalingH,
125                                 final double[] scaled, final Array2DRowRealMatrix nordsieck) {
126            this.referenceTime = referenceTime;
127            this.scalingH      = scalingH;
128            this.scaled        = scaled;
129            this.nordsieck     = nordsieck;
130    
131            // make sure the state and derivatives will depend on the new arrays
132            setInterpolatedTime(getInterpolatedTime());
133    
134        }
135    
136        /** Rescale the instance.
137         * <p>Since the scaled and Nordiseck arrays are shared with the caller,
138         * this method has the side effect of rescaling this arrays in the caller too.</p>
139         * @param scalingH new step size to use in the scaled and nordsieck arrays
140         */
141        public void rescale(final double scalingH) {
142    
143            final double ratio = scalingH / this.scalingH;
144            for (int i = 0; i < scaled.length; ++i) {
145                scaled[i] *= ratio;
146            }
147    
148            final double[][] nData = nordsieck.getDataRef();
149            double power = ratio;
150            for (int i = 0; i < nData.length; ++i) {
151                power *= ratio;
152                final double[] nDataI = nData[i];
153                for (int j = 0; j < nDataI.length; ++j) {
154                    nDataI[j] *= power;
155                }
156            }
157    
158            this.scalingH = scalingH;
159    
160        }
161    
162        /**
163         * Get the state vector variation from current to interpolated state.
164         * <p>This method is aimed at computing y(t<sub>interpolation</sub>)
165         * -y(t<sub>current</sub>) accurately by avoiding the cancellation errors
166         * that would occur if the subtraction were performed explicitly.</p>
167         * <p>The returned vector is a reference to a reused array, so
168         * it should not be modified and it should be copied if it needs
169         * to be preserved across several calls.</p>
170         * @return state vector at time {@link #getInterpolatedTime}
171         * @see #getInterpolatedDerivatives()
172         * @throws DerivativeException if this call induces an automatic
173         * step finalization that throws one
174         */
175        public double[] getInterpolatedStateVariation()
176            throws DerivativeException {
177            // compute and ignore interpolated state
178            // to make sure state variation is computed as a side effect
179            getInterpolatedState();
180            return stateVariation;
181        }
182    
183        /** {@inheritDoc} */
184        @Override
185        protected void computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH) {
186    
187            final double x = interpolatedTime - referenceTime;
188            final double normalizedAbscissa = x / scalingH;
189    
190            Arrays.fill(stateVariation, 0.0);
191            Arrays.fill(interpolatedDerivatives, 0.0);
192    
193            // apply Taylor formula from high order to low order,
194            // for the sake of numerical accuracy
195            final double[][] nData = nordsieck.getDataRef();
196            for (int i = nData.length - 1; i >= 0; --i) {
197                final int order = i + 2;
198                final double[] nDataI = nData[i];
199                final double power = Math.pow(normalizedAbscissa, order);
200                for (int j = 0; j < nDataI.length; ++j) {
201                    final double d = nDataI[j] * power;
202                    stateVariation[j]          += d;
203                    interpolatedDerivatives[j] += order * d;
204                }
205            }
206    
207            for (int j = 0; j < currentState.length; ++j) {
208                stateVariation[j] += scaled[j] * normalizedAbscissa;
209                interpolatedState[j] = currentState[j] + stateVariation[j];
210                interpolatedDerivatives[j] =
211                    (interpolatedDerivatives[j] + scaled[j] * normalizedAbscissa) / x;
212            }
213    
214        }
215    
216        /** {@inheritDoc} */
217        @Override
218        public void writeExternal(final ObjectOutput out)
219            throws IOException {
220    
221            // save the state of the base class
222            writeBaseExternal(out);
223    
224            // save the local attributes
225            out.writeDouble(scalingH);
226            out.writeDouble(referenceTime);
227    
228            final int n = (currentState == null) ? -1 : currentState.length;
229            if (scaled == null) {
230                out.writeBoolean(false);
231            } else {
232                out.writeBoolean(true);
233                for (int j = 0; j < n; ++j) {
234                    out.writeDouble(scaled[j]);
235                }
236            }
237    
238            if (nordsieck == null) {
239                out.writeBoolean(false);
240            } else {
241                out.writeBoolean(true);
242                out.writeObject(nordsieck);
243            }
244    
245            // we don't save state variation, it will be recomputed
246    
247        }
248    
249        /** {@inheritDoc} */
250        @Override
251        public void readExternal(final ObjectInput in)
252            throws IOException, ClassNotFoundException {
253    
254            // read the base class 
255            final double t = readBaseExternal(in);
256    
257            // read the local attributes
258            scalingH      = in.readDouble();
259            referenceTime = in.readDouble();
260    
261            final int n = (currentState == null) ? -1 : currentState.length;
262            final boolean hasScaled = in.readBoolean();
263            if (hasScaled) {
264                scaled = new double[n];
265                for (int j = 0; j < n; ++j) {
266                    scaled[j] = in.readDouble();
267                }
268            } else {
269                scaled = null;
270            }
271    
272            final boolean hasNordsieck = in.readBoolean();
273            if (hasNordsieck) {
274                nordsieck = (Array2DRowRealMatrix) in.readObject();
275            } else {
276                nordsieck = null;
277            }
278    
279            if (hasScaled && hasNordsieck) {
280                // we can now set the interpolated time and state
281                stateVariation = new double[n];
282                setInterpolatedTime(t);
283            } else {
284                stateVariation = null;
285            }
286    
287        }
288    
289    }