/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.ode.nonstiff;
import java.io.IOException;
import java.io.ObjectInput;
import java.io.ObjectOutput;
import org.apache.commons.math3.ode.EquationsMapper;
import org.apache.commons.math3.ode.sampling.AbstractStepInterpolator;
import org.apache.commons.math3.ode.sampling.StepInterpolator;
import org.apache.commons.math3.util.FastMath;
/**
* This class implements an interpolator for the Gragg-Bulirsch-Stoer
* integrator.
*
* <p>This interpolator compute dense output inside the last step
* produced by a Gragg-Bulirsch-Stoer integrator.</p>
*
* <p>
* This implementation is basically a reimplementation in Java of the
* <a
* href="http://www.unige.ch/math/folks/hairer/prog/nonstiff/odex.f">odex</a>
* fortran code by E. Hairer and G. Wanner. The redistribution policy
* for this code is available <a
* href="http://www.unige.ch/~hairer/prog/licence.txt">here</a>, for
* convenience, it is reproduced below.</p>
* </p>
*
* <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
* <tr><td>Copyright (c) 2004, Ernst Hairer</td></tr>
*
* <tr><td>Redistribution and use in source and binary forms, with or
* without modification, are permitted provided that the following
* conditions are met:
* <ul>
* <li>Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.</li>
* <li>Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.</li>
* </ul></td></tr>
*
* <tr><td><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
* CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
* BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR
* CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
* EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
* PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
* PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
* LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
* NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.</strong></td></tr>
* </table>
*
* @see GraggBulirschStoerIntegrator
* @since 1.2
*/
class GraggBulirschStoerStepInterpolator
extends AbstractStepInterpolator {
/** Serializable version identifier. */
private static final long serialVersionUID = 20110928L;
/** Slope at the beginning of the step. */
private double[] y0Dot;
/** State at the end of the step. */
private double[] y1;
/** Slope at the end of the step. */
private double[] y1Dot;
/** Derivatives at the middle of the step.
* element 0 is state at midpoint, element 1 is first derivative ...
*/
private double[][] yMidDots;
/** Interpolation polynomials. */
private double[][] polynomials;
/** Error coefficients for the interpolation. */
private double[] errfac;
/** Degree of the interpolation polynomials. */
private int currentDegree;
/** Simple constructor.
* This constructor should not be used directly, it is only intended
* for the serialization process.
*/
public GraggBulirschStoerStepInterpolator() {
y0Dot = null;
y1 = null;
y1Dot = null;
yMidDots = null;
resetTables(-1);
}
/** Simple constructor.
* @param y reference to the integrator array holding the current state
* @param y0Dot reference to the integrator array holding the slope
* at the beginning of the step
* @param y1 reference to the integrator array holding the state at
* the end of the step
* @param y1Dot reference to the integrator array holding the slope
* at the end of the step
* @param yMidDots reference to the integrator array holding the
* derivatives at the middle point of the step
* @param forward integration direction indicator
* @param primaryMapper equations mapper for the primary equations set
* @param secondaryMappers equations mappers for the secondary equations sets
*/
public GraggBulirschStoerStepInterpolator(final double[] y, final double[] y0Dot,
final double[] y1, final double[] y1Dot,
final double[][] yMidDots,
final boolean forward,
final EquationsMapper primaryMapper,
final EquationsMapper[] secondaryMappers) {
super(y, forward, primaryMapper, secondaryMappers);
this.y0Dot = y0Dot;
this.y1 = y1;
this.y1Dot = y1Dot;
this.yMidDots = yMidDots;
resetTables(yMidDots.length + 4);
}
/** Copy constructor.
* @param interpolator interpolator to copy from. The copy is a deep
* copy: its arrays are separated from the original arrays of the
* instance
*/
public GraggBulirschStoerStepInterpolator
(final GraggBulirschStoerStepInterpolator interpolator) {
super(interpolator);
final int dimension = currentState.length;
// the interpolator has been finalized,
// the following arrays are not needed anymore
y0Dot = null;
y1 = null;
y1Dot = null;
yMidDots = null;
// copy the interpolation polynomials (up to the current degree only)
if (interpolator.polynomials == null) {
polynomials = null;
currentDegree = -1;
} else {
resetTables(interpolator.currentDegree);
for (int i = 0; i < polynomials.length; ++i) {
polynomials[i] = new double[dimension];
System.arraycopy(interpolator.polynomials[i], 0,
polynomials[i], 0, dimension);
}
currentDegree = interpolator.currentDegree;
}
}
/** Reallocate the internal tables.
* Reallocate the internal tables in order to be able to handle
* interpolation polynomials up to the given degree
* @param maxDegree maximal degree to handle
*/
private void resetTables(final int maxDegree) {
if (maxDegree < 0) {
polynomials = null;
errfac = null;
currentDegree = -1;
} else {
final double[][] newPols = new double[maxDegree + 1][];
if (polynomials != null) {
System.arraycopy(polynomials, 0, newPols, 0, polynomials.length);
for (int i = polynomials.length; i < newPols.length; ++i) {
newPols[i] = new double[currentState.length];
}
} else {
for (int i = 0; i < newPols.length; ++i) {
newPols[i] = new double[currentState.length];
}
}
polynomials = newPols;
// initialize the error factors array for interpolation
if (maxDegree <= 4) {
errfac = null;
} else {
errfac = new double[maxDegree - 4];
for (int i = 0; i < errfac.length; ++i) {
final int ip5 = i + 5;
errfac[i] = 1.0 / (ip5 * ip5);
final double e = 0.5 * FastMath.sqrt (((double) (i + 1)) / ip5);
for (int j = 0; j <= i; ++j) {
errfac[i] *= e / (j + 1);
}
}
}
currentDegree = 0;
}
}
/** {@inheritDoc} */
@Override
protected StepInterpolator doCopy() {
return new GraggBulirschStoerStepInterpolator(this);
}
/** Compute the interpolation coefficients for dense output.
* @param mu degree of the interpolation polynomial
* @param h current step
*/
public void computeCoefficients(final int mu, final double h) {
if ((polynomials == null) || (polynomials.length <= (mu + 4))) {
resetTables(mu + 4);
}
currentDegree = mu + 4;
for (int i = 0; i < currentState.length; ++i) {
final double yp0 = h * y0Dot[i];
final double yp1 = h * y1Dot[i];
final double ydiff = y1[i] - currentState[i];
final double aspl = ydiff - yp1;
final double bspl = yp0 - ydiff;
polynomials[0][i] = currentState[i];
polynomials[1][i] = ydiff;
polynomials[2][i] = aspl;
polynomials[3][i] = bspl;
if (mu < 0) {
return;
}
// compute the remaining coefficients
final double ph0 = 0.5 * (currentState[i] + y1[i]) + 0.125 * (aspl + bspl);
polynomials[4][i] = 16 * (yMidDots[0][i] - ph0);
if (mu > 0) {
final double ph1 = ydiff + 0.25 * (aspl - bspl);
polynomials[5][i] = 16 * (yMidDots[1][i] - ph1);
if (mu > 1) {
final double ph2 = yp1 - yp0;
polynomials[6][i] = 16 * (yMidDots[2][i] - ph2 + polynomials[4][i]);
if (mu > 2) {
final double ph3 = 6 * (bspl - aspl);
polynomials[7][i] = 16 * (yMidDots[3][i] - ph3 + 3 * polynomials[5][i]);
for (int j = 4; j <= mu; ++j) {
final double fac1 = 0.5 * j * (j - 1);
final double fac2 = 2 * fac1 * (j - 2) * (j - 3);
polynomials[j+4][i] =
16 * (yMidDots[j][i] + fac1 * polynomials[j+2][i] - fac2 * polynomials[j][i]);
}
}
}
}
}
}
/** Estimate interpolation error.
* @param scale scaling array
* @return estimate of the interpolation error
*/
public double estimateError(final double[] scale) {
double error = 0;
if (currentDegree >= 5) {
for (int i = 0; i < scale.length; ++i) {
final double e = polynomials[currentDegree][i] / scale[i];
error += e * e;
}
error = FastMath.sqrt(error / scale.length) * errfac[currentDegree - 5];
}
return error;
}
/** {@inheritDoc} */
@Override
protected void computeInterpolatedStateAndDerivatives(final double theta,
final double oneMinusThetaH) {
final int dimension = currentState.length;
final double oneMinusTheta = 1.0 - theta;
final double theta05 = theta - 0.5;
final double tOmT = theta * oneMinusTheta;
final double t4 = tOmT * tOmT;
final double t4Dot = 2 * tOmT * (1 - 2 * theta);
final double dot1 = 1.0 / h;
final double dot2 = theta * (2 - 3 * theta) / h;
final double dot3 = ((3 * theta - 4) * theta + 1) / h;
for (int i = 0; i < dimension; ++i) {
final double p0 = polynomials[0][i];
final double p1 = polynomials[1][i];
final double p2 = polynomials[2][i];
final double p3 = polynomials[3][i];
interpolatedState[i] = p0 + theta * (p1 + oneMinusTheta * (p2 * theta + p3 * oneMinusTheta));
interpolatedDerivatives[i] = dot1 * p1 + dot2 * p2 + dot3 * p3;
if (currentDegree > 3) {
double cDot = 0;
double c = polynomials[currentDegree][i];
for (int j = currentDegree - 1; j > 3; --j) {
final double d = 1.0 / (j - 3);
cDot = d * (theta05 * cDot + c);
c = polynomials[j][i] + c * d * theta05;
}
interpolatedState[i] += t4 * c;
interpolatedDerivatives[i] += (t4 * cDot + t4Dot * c) / h;
}
}
if (h == 0) {
// in this degenerated case, the previous computation leads to NaN for derivatives
// we fix this by using the derivatives at midpoint
System.arraycopy(yMidDots[1], 0, interpolatedDerivatives, 0, dimension);
}
}
/** {@inheritDoc} */
@Override
public void writeExternal(final ObjectOutput out)
throws IOException {
final int dimension = (currentState == null) ? -1 : currentState.length;
// save the state of the base class
writeBaseExternal(out);
// save the local attributes (but not the temporary vectors)
out.writeInt(currentDegree);
for (int k = 0; k <= currentDegree; ++k) {
for (int l = 0; l < dimension; ++l) {
out.writeDouble(polynomials[k][l]);
}
}
}
/** {@inheritDoc} */
@Override
public void readExternal(final ObjectInput in)
throws IOException, ClassNotFoundException {
// read the base class
final double t = readBaseExternal(in);
final int dimension = (currentState == null) ? -1 : currentState.length;
// read the local attributes
final int degree = in.readInt();
resetTables(degree);
currentDegree = degree;
for (int k = 0; k <= currentDegree; ++k) {
for (int l = 0; l < dimension; ++l) {
polynomials[k][l] = in.readDouble();
}
}
// we can now set the interpolated time and state
setInterpolatedTime(t);
}
}