View Javadoc

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  
18  package org.apache.commons.math.ode.sampling;
19  
20  import java.io.IOException;
21  import java.io.ObjectInput;
22  import java.io.ObjectOutput;
23  import java.util.Arrays;
24  
25  import org.apache.commons.math.linear.Array2DRowRealMatrix;
26  import org.apache.commons.math.ode.DerivativeException;
27  
28  /**
29   * This class implements an interpolator for integrators using Nordsieck representation.
30   *
31   * <p>This interpolator computes dense output around the current point.
32   * The interpolation equation is based on Taylor series formulas.
33   *
34   * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator
35   * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator
36   * @version $Revision: 791723 $ $Date: 2009-07-07 03:07:23 -0400 (Tue, 07 Jul 2009) $
37   * @since 2.0
38   */
39  
40  public class NordsieckStepInterpolator extends AbstractStepInterpolator {
41  
42      /** Serializable version identifier */
43      private static final long serialVersionUID = -7179861704951334960L;
44  
45      /** Step size used in the first scaled derivative and Nordsieck vector. */
46      private double scalingH;
47  
48      /** Reference time for all arrays.
49       * <p>Sometimes, the reference time is the same as previousTime,
50       * sometimes it is the same as currentTime, so we use a separate
51       * field to avoid any confusion.
52       * </p>
53       */
54      private double referenceTime;
55  
56      /** First scaled derivative. */
57      private double[] scaled;
58  
59      /** Nordsieck vector. */
60      private Array2DRowRealMatrix nordsieck;
61  
62      /** State variation. */
63      protected double[] stateVariation;
64  
65      /** Simple constructor.
66       * This constructor builds an instance that is not usable yet, the
67       * {@link AbstractStepInterpolator#reinitialize} method should be called
68       * before using the instance in order to initialize the internal arrays. This
69       * constructor is used only in order to delay the initialization in
70       * some cases.
71       */
72      public NordsieckStepInterpolator() {
73      }
74  
75      /** Copy constructor.
76       * @param interpolator interpolator to copy from. The copy is a deep
77       * copy: its arrays are separated from the original arrays of the
78       * instance
79       */
80      public NordsieckStepInterpolator(final NordsieckStepInterpolator interpolator) {
81          super(interpolator);
82          scalingH      = interpolator.scalingH;
83          referenceTime = interpolator.referenceTime;
84          if (interpolator.scaled != null) {
85              scaled = interpolator.scaled.clone();
86          }
87          if (interpolator.nordsieck != null) {
88              nordsieck = new Array2DRowRealMatrix(interpolator.nordsieck.getDataRef(), true);
89          }
90          if (interpolator.stateVariation != null) {
91              stateVariation = interpolator.stateVariation.clone();
92          }
93      }
94  
95      /** {@inheritDoc} */
96      @Override
97      protected StepInterpolator doCopy() {
98          return new NordsieckStepInterpolator(this);
99      }
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 }