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.math3.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.math3.exception.MaxCountExceededException;
026    import org.apache.commons.math3.linear.Array2DRowRealMatrix;
027    import org.apache.commons.math3.ode.EquationsMapper;
028    import 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 1416643 2012-12-03 19:37:14Z tn $
039     * @since 2.0
040     */
041    
042    public 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 Nordiseck 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    }