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.nonstiff;
19  
20  import java.io.IOException;
21  import java.io.ObjectInput;
22  import java.io.ObjectOutput;
23  
24  import org.apache.commons.math.ode.DerivativeException;
25  import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
26  import org.apache.commons.math.ode.sampling.StepInterpolator;
27  
28  /**
29   * This class implements an interpolator for the Gragg-Bulirsch-Stoer
30   * integrator.
31   *
32   * <p>This interpolator compute dense output inside the last step
33   * produced by a Gragg-Bulirsch-Stoer integrator.</p>
34   *
35   * <p>
36   * This implementation is basically a reimplementation in Java of the
37   * <a
38   * href="http://www.unige.ch/math/folks/hairer/prog/nonstiff/odex.f">odex</a>
39   * fortran code by E. Hairer and G. Wanner. The redistribution policy
40   * for this code is available <a
41   * href="http://www.unige.ch/~hairer/prog/licence.txt">here</a>, for
42   * convenience, it is reproduced below.</p>
43   * </p>
44   *
45   * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
46   * <tr><td>Copyright (c) 2004, Ernst Hairer</td></tr>
47   *
48   * <tr><td>Redistribution and use in source and binary forms, with or
49   * without modification, are permitted provided that the following
50   * conditions are met:
51   * <ul>
52   *  <li>Redistributions of source code must retain the above copyright
53   *      notice, this list of conditions and the following disclaimer.</li>
54   *  <li>Redistributions in binary form must reproduce the above copyright
55   *      notice, this list of conditions and the following disclaimer in the
56   *      documentation and/or other materials provided with the distribution.</li>
57   * </ul></td></tr>
58   *
59   * <tr><td><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
60   * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
61   * BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
62   * FOR A  PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR
63   * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
64   * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
65   * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
66   * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
67   * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
68   * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
69   * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.</strong></td></tr>
70   * </table>
71   *
72   * @see GraggBulirschStoerIntegrator
73   * @version $Revision: 782432 $ $Date: 2009-06-07 15:08:26 -0400 (Sun, 07 Jun 2009) $
74   * @author E. Hairer and G. Wanner (fortran version)
75   * @since 1.2
76   */
77  
78  class GraggBulirschStoerStepInterpolator
79    extends AbstractStepInterpolator {
80  
81    /** Slope at the beginning of the step. */
82    private double[] y0Dot;
83  
84    /** State at the end of the step. */
85    private double[] y1;
86  
87    /** Slope at the end of the step. */
88    private double[] y1Dot;
89  
90    /** Derivatives at the middle of the step.
91     * element 0 is state at midpoint, element 1 is first derivative ...
92     */
93    private double[][] yMidDots;
94  
95    /** Interpolation polynoms. */
96    private double[][] polynoms;
97  
98    /** Error coefficients for the interpolation. */
99    private double[] errfac;
100 
101   /** Degree of the interpolation polynoms. */
102   private int currentDegree;
103 
104   /** Reallocate the internal tables.
105    * Reallocate the internal tables in order to be able to handle
106    * interpolation polynoms up to the given degree
107    * @param maxDegree maximal degree to handle
108    */
109   private void resetTables(final int maxDegree) {
110 
111     if (maxDegree < 0) {
112       polynoms      = null;
113       errfac        = null;
114       currentDegree = -1;
115     } else {
116 
117       final double[][] newPols = new double[maxDegree + 1][];
118       if (polynoms != null) {
119         System.arraycopy(polynoms, 0, newPols, 0, polynoms.length);
120         for (int i = polynoms.length; i < newPols.length; ++i) {
121           newPols[i] = new double[currentState.length];
122         }
123       } else {
124         for (int i = 0; i < newPols.length; ++i) {
125           newPols[i] = new double[currentState.length];
126         }
127       }
128       polynoms = newPols;
129 
130       // initialize the error factors array for interpolation
131       if (maxDegree <= 4) {
132         errfac = null;
133       } else {
134         errfac = new double[maxDegree - 4];
135         for (int i = 0; i < errfac.length; ++i) {
136           final int ip5 = i + 5;
137           errfac[i] = 1.0 / (ip5 * ip5);
138           final double e = 0.5 * Math.sqrt (((double) (i + 1)) / ip5);
139           for (int j = 0; j <= i; ++j) {
140             errfac[i] *= e / (j + 1);
141           }
142         }
143       }
144 
145       currentDegree = 0;
146 
147     }
148 
149   }
150 
151   /** Simple constructor.
152    * This constructor should not be used directly, it is only intended
153    * for the serialization process.
154    */
155   public GraggBulirschStoerStepInterpolator() {
156     y0Dot    = null;
157     y1       = null;
158     y1Dot    = null;
159     yMidDots = null;
160     resetTables(-1);
161   }
162 
163   /** Simple constructor.
164    * @param y reference to the integrator array holding the current state
165    * @param y0Dot reference to the integrator array holding the slope
166    * at the beginning of the step
167    * @param y1 reference to the integrator array holding the state at
168    * the end of the step
169    * @param y1Dot reference to the integrator array holding the slope
170    * at theend of the step
171    * @param yMidDots reference to the integrator array holding the
172    * derivatives at the middle point of the step
173    * @param forward integration direction indicator
174    */
175   public GraggBulirschStoerStepInterpolator(final double[] y, final double[] y0Dot,
176                                             final double[] y1, final double[] y1Dot,
177                                             final double[][] yMidDots,
178                                             final boolean forward) {
179 
180     super(y, forward);
181     this.y0Dot    = y0Dot;
182     this.y1       = y1;
183     this.y1Dot    = y1Dot;
184     this.yMidDots = yMidDots;
185 
186     resetTables(yMidDots.length + 4);
187 
188   }
189 
190   /** Copy constructor.
191    * @param interpolator interpolator to copy from. The copy is a deep
192    * copy: its arrays are separated from the original arrays of the
193    * instance
194    */
195   public GraggBulirschStoerStepInterpolator
196     (final GraggBulirschStoerStepInterpolator interpolator) {
197 
198     super(interpolator);
199 
200     final int dimension = currentState.length;
201 
202     // the interpolator has been finalized,
203     // the following arrays are not needed anymore
204     y0Dot    = null;
205     y1       = null;
206     y1Dot    = null;
207     yMidDots = null;
208 
209     // copy the interpolation polynoms (up to the current degree only)
210     if (interpolator.polynoms == null) {
211       polynoms = null;
212       currentDegree = -1;
213     } else {
214       resetTables(interpolator.currentDegree);
215       for (int i = 0; i < polynoms.length; ++i) {
216         polynoms[i] = new double[dimension];
217         System.arraycopy(interpolator.polynoms[i], 0,
218                          polynoms[i], 0, dimension);
219       }
220       currentDegree = interpolator.currentDegree;
221     }
222 
223   }
224 
225   /** {@inheritDoc} */
226   @Override
227   protected StepInterpolator doCopy() {
228     return new GraggBulirschStoerStepInterpolator(this);
229   }
230 
231 
232   /** Compute the interpolation coefficients for dense output.
233    * @param mu degree of the interpolation polynomial
234    * @param h current step
235    */
236   public void computeCoefficients(final int mu, final double h) {
237 
238     if ((polynoms == null) || (polynoms.length <= (mu + 4))) {
239       resetTables(mu + 4);
240     }
241 
242     currentDegree = mu + 4;
243 
244     for (int i = 0; i < currentState.length; ++i) {
245 
246       final double yp0   = h * y0Dot[i];
247       final double yp1   = h * y1Dot[i];
248       final double ydiff = y1[i] - currentState[i];
249       final double aspl  = ydiff - yp1;
250       final double bspl  = yp0 - ydiff;
251 
252       polynoms[0][i] = currentState[i];
253       polynoms[1][i] = ydiff;
254       polynoms[2][i] = aspl;
255       polynoms[3][i] = bspl;
256 
257       if (mu < 0) {
258         return;
259       }
260 
261       // compute the remaining coefficients
262       final double ph0 = 0.5 * (currentState[i] + y1[i]) + 0.125 * (aspl + bspl);
263       polynoms[4][i] = 16 * (yMidDots[0][i] - ph0);
264 
265       if (mu > 0) {
266         final double ph1 = ydiff + 0.25 * (aspl - bspl);
267         polynoms[5][i] = 16 * (yMidDots[1][i] - ph1);
268 
269         if (mu > 1) {
270           final double ph2 = yp1 - yp0;
271           polynoms[6][i] = 16 * (yMidDots[2][i] - ph2 + polynoms[4][i]);
272 
273           if (mu > 2) {
274             final double ph3 = 6 * (bspl - aspl);
275             polynoms[7][i] = 16 * (yMidDots[3][i] - ph3 + 3 * polynoms[5][i]);
276 
277             for (int j = 4; j <= mu; ++j) {
278               final double fac1 = 0.5 * j * (j - 1);
279               final double fac2 = 2 * fac1 * (j - 2) * (j - 3);
280               polynoms[j+4][i] =
281                   16 * (yMidDots[j][i] + fac1 * polynoms[j+2][i] - fac2 * polynoms[j][i]);
282             }
283 
284           }
285         }
286       }
287     }
288 
289   }
290 
291   /** Estimate interpolation error.
292    * @param scale scaling array
293    * @return estimate of the interpolation error
294    */
295   public double estimateError(final double[] scale) {
296     double error = 0;
297     if (currentDegree >= 5) {
298       for (int i = 0; i < currentState.length; ++i) {
299         final double e = polynoms[currentDegree][i] / scale[i];
300         error += e * e;
301       }
302       error = Math.sqrt(error / currentState.length) * errfac[currentDegree-5];
303     }
304     return error;
305   }
306 
307   /** {@inheritDoc} */
308   @Override
309   protected void computeInterpolatedStateAndDerivatives(final double theta,
310                                           final double oneMinusThetaH)
311     throws DerivativeException {
312 
313     final int dimension = currentState.length;
314 
315     final double oneMinusTheta = 1.0 - theta;
316     final double theta05       = theta - 0.5;
317     final double tOmT          = theta * oneMinusTheta;
318     final double t4            = tOmT * tOmT;
319     final double t4Dot         = 2 * tOmT * (1 - 2 * theta);
320     final double dot1          = 1.0 / h;
321     final double dot2          = theta * (2 - 3 * theta) / h;
322     final double dot3          = ((3 * theta - 4) * theta + 1) / h;
323 
324     for (int i = 0; i < dimension; ++i) {
325 
326         final double p0 = polynoms[0][i];
327         final double p1 = polynoms[1][i];
328         final double p2 = polynoms[2][i];
329         final double p3 = polynoms[3][i];
330         interpolatedState[i] = p0 + theta * (p1 + oneMinusTheta * (p2 * theta + p3 * oneMinusTheta));
331         interpolatedDerivatives[i] = dot1 * p1 + dot2 * p2 + dot3 * p3;
332 
333         if (currentDegree > 3) {
334             double cDot = 0;
335             double c = polynoms[currentDegree][i];
336             for (int j = currentDegree - 1; j > 3; --j) {
337                 final double d = 1.0 / (j - 3);
338                 cDot = d * (theta05 * cDot + c);
339                 c = polynoms[j][i] + c * d * theta05;
340             }
341             interpolatedState[i]       += t4 * c;
342             interpolatedDerivatives[i] += (t4 * cDot + t4Dot * c) / h;
343         }
344 
345     }
346 
347     if (h == 0) {
348         // in this degenerated case, the previous computation leads to NaN for derivatives
349         // we fix this by using the derivatives at midpoint
350         System.arraycopy(yMidDots[1], 0, interpolatedDerivatives, 0, dimension);
351     }
352 
353   }
354     
355   /** {@inheritDoc} */
356   @Override
357   public void writeExternal(final ObjectOutput out)
358     throws IOException {
359 
360     final int dimension = (currentState == null) ? -1 : currentState.length;
361 
362     // save the state of the base class
363     writeBaseExternal(out);
364 
365     // save the local attributes (but not the temporary vectors)
366     out.writeInt(currentDegree);
367     for (int k = 0; k <= currentDegree; ++k) {
368       for (int l = 0; l < dimension; ++l) {
369         out.writeDouble(polynoms[k][l]);
370       }
371     }
372 
373   }
374 
375   /** {@inheritDoc} */
376   @Override
377   public void readExternal(final ObjectInput in)
378     throws IOException {
379 
380     // read the base class 
381     final double t = readBaseExternal(in);
382     final int dimension = (currentState == null) ? -1 : currentState.length;
383 
384     // read the local attributes
385     final int degree = in.readInt();
386     resetTables(degree);
387     currentDegree = degree;
388 
389     for (int k = 0; k <= currentDegree; ++k) {
390       for (int l = 0; l < dimension; ++l) {
391         polynoms[k][l] = in.readDouble();
392       }
393     }
394 
395     // we can now set the interpolated time and state
396     setInterpolatedTime(t);
397 
398   }
399 
400   /** Serializable version identifier */
401   private static final long serialVersionUID = 7320613236731409847L;
402 
403 }