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
18package org.apache.commons.math.ode.nonstiff;
19
20import java.io.IOException;
21import java.io.ObjectInput;
22import java.io.ObjectOutput;
23
24import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
25import org.apache.commons.math.ode.sampling.StepInterpolator;
26import org.apache.commons.math.util.FastMath;
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: 1061507 $ $Date: 2011-01-20 21:55:00 +0100 (jeu. 20 janv. 2011) $
74 * @since 1.2
75 */
76
77class GraggBulirschStoerStepInterpolator
78  extends AbstractStepInterpolator {
79
80    /** Serializable version identifier. */
81    private static final long serialVersionUID = 7320613236731409847L;
82
83    /** Slope at the beginning of the step. */
84    private double[] y0Dot;
85
86    /** State at the end of the step. */
87    private double[] y1;
88
89    /** Slope at the end of the step. */
90    private double[] y1Dot;
91
92    /** Derivatives at the middle of the step.
93     * element 0 is state at midpoint, element 1 is first derivative ...
94     */
95    private double[][] yMidDots;
96
97    /** Interpolation polynoms. */
98    private double[][] polynoms;
99
100    /** Error coefficients for the interpolation. */
101    private double[] errfac;
102
103    /** Degree of the interpolation polynoms. */
104    private int currentDegree;
105
106  /** Simple constructor.
107   * This constructor should not be used directly, it is only intended
108   * for the serialization process.
109   */
110  public GraggBulirschStoerStepInterpolator() {
111    y0Dot    = null;
112    y1       = null;
113    y1Dot    = null;
114    yMidDots = null;
115    resetTables(-1);
116  }
117
118  /** Simple constructor.
119   * @param y reference to the integrator array holding the current state
120   * @param y0Dot reference to the integrator array holding the slope
121   * at the beginning of the step
122   * @param y1 reference to the integrator array holding the state at
123   * the end of the step
124   * @param y1Dot reference to the integrator array holding the slope
125   * at the end of the step
126   * @param yMidDots reference to the integrator array holding the
127   * derivatives at the middle point of the step
128   * @param forward integration direction indicator
129   */
130  public GraggBulirschStoerStepInterpolator(final double[] y, final double[] y0Dot,
131                                            final double[] y1, final double[] y1Dot,
132                                            final double[][] yMidDots,
133                                            final boolean forward) {
134
135    super(y, forward);
136    this.y0Dot    = y0Dot;
137    this.y1       = y1;
138    this.y1Dot    = y1Dot;
139    this.yMidDots = yMidDots;
140
141    resetTables(yMidDots.length + 4);
142
143  }
144
145  /** Copy constructor.
146   * @param interpolator interpolator to copy from. The copy is a deep
147   * copy: its arrays are separated from the original arrays of the
148   * instance
149   */
150  public GraggBulirschStoerStepInterpolator
151    (final GraggBulirschStoerStepInterpolator interpolator) {
152
153    super(interpolator);
154
155    final int dimension = currentState.length;
156
157    // the interpolator has been finalized,
158    // the following arrays are not needed anymore
159    y0Dot    = null;
160    y1       = null;
161    y1Dot    = null;
162    yMidDots = null;
163
164    // copy the interpolation polynoms (up to the current degree only)
165    if (interpolator.polynoms == null) {
166      polynoms = null;
167      currentDegree = -1;
168    } else {
169      resetTables(interpolator.currentDegree);
170      for (int i = 0; i < polynoms.length; ++i) {
171        polynoms[i] = new double[dimension];
172        System.arraycopy(interpolator.polynoms[i], 0,
173                         polynoms[i], 0, dimension);
174      }
175      currentDegree = interpolator.currentDegree;
176    }
177
178  }
179
180  /** Reallocate the internal tables.
181   * Reallocate the internal tables in order to be able to handle
182   * interpolation polynoms up to the given degree
183   * @param maxDegree maximal degree to handle
184   */
185  private void resetTables(final int maxDegree) {
186
187    if (maxDegree < 0) {
188      polynoms      = null;
189      errfac        = null;
190      currentDegree = -1;
191    } else {
192
193      final double[][] newPols = new double[maxDegree + 1][];
194      if (polynoms != null) {
195        System.arraycopy(polynoms, 0, newPols, 0, polynoms.length);
196        for (int i = polynoms.length; i < newPols.length; ++i) {
197          newPols[i] = new double[currentState.length];
198        }
199      } else {
200        for (int i = 0; i < newPols.length; ++i) {
201          newPols[i] = new double[currentState.length];
202        }
203      }
204      polynoms = newPols;
205
206      // initialize the error factors array for interpolation
207      if (maxDegree <= 4) {
208        errfac = null;
209      } else {
210        errfac = new double[maxDegree - 4];
211        for (int i = 0; i < errfac.length; ++i) {
212          final int ip5 = i + 5;
213          errfac[i] = 1.0 / (ip5 * ip5);
214          final double e = 0.5 * FastMath.sqrt (((double) (i + 1)) / ip5);
215          for (int j = 0; j <= i; ++j) {
216            errfac[i] *= e / (j + 1);
217          }
218        }
219      }
220
221      currentDegree = 0;
222
223    }
224
225  }
226
227  /** {@inheritDoc} */
228  @Override
229  protected StepInterpolator doCopy() {
230    return new GraggBulirschStoerStepInterpolator(this);
231  }
232
233
234  /** Compute the interpolation coefficients for dense output.
235   * @param mu degree of the interpolation polynomial
236   * @param h current step
237   */
238  public void computeCoefficients(final int mu, final double h) {
239
240    if ((polynoms == null) || (polynoms.length <= (mu + 4))) {
241      resetTables(mu + 4);
242    }
243
244    currentDegree = mu + 4;
245
246    for (int i = 0; i < currentState.length; ++i) {
247
248      final double yp0   = h * y0Dot[i];
249      final double yp1   = h * y1Dot[i];
250      final double ydiff = y1[i] - currentState[i];
251      final double aspl  = ydiff - yp1;
252      final double bspl  = yp0 - ydiff;
253
254      polynoms[0][i] = currentState[i];
255      polynoms[1][i] = ydiff;
256      polynoms[2][i] = aspl;
257      polynoms[3][i] = bspl;
258
259      if (mu < 0) {
260        return;
261      }
262
263      // compute the remaining coefficients
264      final double ph0 = 0.5 * (currentState[i] + y1[i]) + 0.125 * (aspl + bspl);
265      polynoms[4][i] = 16 * (yMidDots[0][i] - ph0);
266
267      if (mu > 0) {
268        final double ph1 = ydiff + 0.25 * (aspl - bspl);
269        polynoms[5][i] = 16 * (yMidDots[1][i] - ph1);
270
271        if (mu > 1) {
272          final double ph2 = yp1 - yp0;
273          polynoms[6][i] = 16 * (yMidDots[2][i] - ph2 + polynoms[4][i]);
274
275          if (mu > 2) {
276            final double ph3 = 6 * (bspl - aspl);
277            polynoms[7][i] = 16 * (yMidDots[3][i] - ph3 + 3 * polynoms[5][i]);
278
279            for (int j = 4; j <= mu; ++j) {
280              final double fac1 = 0.5 * j * (j - 1);
281              final double fac2 = 2 * fac1 * (j - 2) * (j - 3);
282              polynoms[j+4][i] =
283                  16 * (yMidDots[j][i] + fac1 * polynoms[j+2][i] - fac2 * polynoms[j][i]);
284            }
285
286          }
287        }
288      }
289    }
290
291  }
292
293  /** Estimate interpolation error.
294   * @param scale scaling array
295   * @return estimate of the interpolation error
296   */
297  public double estimateError(final double[] scale) {
298    double error = 0;
299    if (currentDegree >= 5) {
300      for (int i = 0; i < scale.length; ++i) {
301        final double e = polynoms[currentDegree][i] / scale[i];
302        error += e * e;
303      }
304      error = FastMath.sqrt(error / scale.length) * errfac[currentDegree - 5];
305    }
306    return error;
307  }
308
309  /** {@inheritDoc} */
310  @Override
311  protected void computeInterpolatedStateAndDerivatives(final double theta,
312                                          final double oneMinusThetaH) {
313
314    final int dimension = currentState.length;
315
316    final double oneMinusTheta = 1.0 - theta;
317    final double theta05       = theta - 0.5;
318    final double tOmT          = theta * oneMinusTheta;
319    final double t4            = tOmT * tOmT;
320    final double t4Dot         = 2 * tOmT * (1 - 2 * theta);
321    final double dot1          = 1.0 / h;
322    final double dot2          = theta * (2 - 3 * theta) / h;
323    final double dot3          = ((3 * theta - 4) * theta + 1) / h;
324
325    for (int i = 0; i < dimension; ++i) {
326
327        final double p0 = polynoms[0][i];
328        final double p1 = polynoms[1][i];
329        final double p2 = polynoms[2][i];
330        final double p3 = polynoms[3][i];
331        interpolatedState[i] = p0 + theta * (p1 + oneMinusTheta * (p2 * theta + p3 * oneMinusTheta));
332        interpolatedDerivatives[i] = dot1 * p1 + dot2 * p2 + dot3 * p3;
333
334        if (currentDegree > 3) {
335            double cDot = 0;
336            double c = polynoms[currentDegree][i];
337            for (int j = currentDegree - 1; j > 3; --j) {
338                final double d = 1.0 / (j - 3);
339                cDot = d * (theta05 * cDot + c);
340                c = polynoms[j][i] + c * d * theta05;
341            }
342            interpolatedState[i]       += t4 * c;
343            interpolatedDerivatives[i] += (t4 * cDot + t4Dot * c) / h;
344        }
345
346    }
347
348    if (h == 0) {
349        // in this degenerated case, the previous computation leads to NaN for derivatives
350        // we fix this by using the derivatives at midpoint
351        System.arraycopy(yMidDots[1], 0, interpolatedDerivatives, 0, dimension);
352    }
353
354  }
355
356  /** {@inheritDoc} */
357  @Override
358  public void writeExternal(final ObjectOutput out)
359    throws IOException {
360
361    final int dimension = (currentState == null) ? -1 : currentState.length;
362
363    // save the state of the base class
364    writeBaseExternal(out);
365
366    // save the local attributes (but not the temporary vectors)
367    out.writeInt(currentDegree);
368    for (int k = 0; k <= currentDegree; ++k) {
369      for (int l = 0; l < dimension; ++l) {
370        out.writeDouble(polynoms[k][l]);
371      }
372    }
373
374  }
375
376  /** {@inheritDoc} */
377  @Override
378  public void readExternal(final ObjectInput in)
379    throws IOException {
380
381    // read the base class
382    final double t = readBaseExternal(in);
383    final int dimension = (currentState == null) ? -1 : currentState.length;
384
385    // read the local attributes
386    final int degree = in.readInt();
387    resetTables(degree);
388    currentDegree = degree;
389
390    for (int k = 0; k <= currentDegree; ++k) {
391      for (int l = 0; l < dimension; ++l) {
392        polynoms[k][l] = in.readDouble();
393      }
394    }
395
396    // we can now set the interpolated time and state
397    setInterpolatedTime(t);
398
399  }
400
401}
402