alvinalexander.com | career | drupal | java | mac | mysql | perl | scala | uml | unix  
* * <tr> * * <tr> * </table> * * @see GraggBulirschStoerIntegrator * @version $Revision: 919479 $ $Date: 2010-03-05 11:35:56 -0500 (Fri, 05 Mar 2010) $ * @since 1.2 */ class GraggBulirschStoerStepInterpolator extends AbstractStepInterpolator { /** Serializable version identifier. */ private static final long serialVersionUID = 7320613236731409847L; /** 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 polynoms. */ private double[][] polynoms; /** Error coefficients for the interpolation. */ private double[] errfac; /** Degree of the interpolation polynoms. */ 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 */ public GraggBulirschStoerStepInterpolator(final double[] y, final double[] y0Dot, final double[] y1, final double[] y1Dot, final double[][] yMidDots, final boolean forward) { super(y, forward); 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 polynoms (up to the current degree only) if (interpolator.polynoms == null) { polynoms = null; currentDegree = -1; } else { resetTables(interpolator.currentDegree); for (int i = 0; i < polynoms.length; ++i) { polynoms[i] = new double[dimension]; System.arraycopy(interpolator.polynoms[i], 0, polynoms[i], 0, dimension); } currentDegree = interpolator.currentDegree; } } /** Reallocate the internal tables. * Reallocate the internal tables in order to be able to handle * interpolation polynoms up to the given degree * @param maxDegree maximal degree to handle */ private void resetTables(final int maxDegree) { if (maxDegree < 0) { polynoms = null; errfac = null; currentDegree = -1; } else { final double[][] newPols = new double[maxDegree + 1][]; if (polynoms != null) { System.arraycopy(polynoms, 0, newPols, 0, polynoms.length); for (int i = polynoms.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]; } } polynoms = 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 * Math.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 ((polynoms == null) || (polynoms.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; polynoms[0][i] = currentState[i]; polynoms[1][i] = ydiff; polynoms[2][i] = aspl; polynoms[3][i] = bspl; if (mu < 0) { return; } // compute the remaining coefficients final double ph0 = 0.5 * (currentState[i] + y1[i]) + 0.125 * (aspl + bspl); polynoms[4][i] = 16 * (yMidDots[0][i] - ph0); if (mu > 0) { final double ph1 = ydiff + 0.25 * (aspl - bspl); polynoms[5][i] = 16 * (yMidDots[1][i] - ph1); if (mu > 1) { final double ph2 = yp1 - yp0; polynoms[6][i] = 16 * (yMidDots[2][i] - ph2 + polynoms[4][i]); if (mu > 2) { final double ph3 = 6 * (bspl - aspl); polynoms[7][i] = 16 * (yMidDots[3][i] - ph3 + 3 * polynoms[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); polynoms[j+4][i] = 16 * (yMidDots[j][i] + fac1 * polynoms[j+2][i] - fac2 * polynoms[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 < currentState.length; ++i) { final double e = polynoms[currentDegree][i] / scale[i]; error += e * e; } error = Math.sqrt(error / currentState.length) * errfac[currentDegree-5]; } return error; } /** {@inheritDoc} */ @Override protected void computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH) throws DerivativeException { 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 = polynoms[0][i]; final double p1 = polynoms[1][i]; final double p2 = polynoms[2][i]; final double p3 = polynoms[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 = polynoms[currentDegree][i]; for (int j = currentDegree - 1; j > 3; --j) { final double d = 1.0 / (j - 3); cDot = d * (theta05 * cDot + c); c = polynoms[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(polynoms[k][l]); } } } /** {@inheritDoc} */ @Override public void readExternal(final ObjectInput in) throws IOException { // 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) { polynoms[k][l] = in.readDouble(); } } // we can now set the interpolated time and state setInterpolatedTime(t); } }

Other Commons Math examples (source code examples)

Here is a short list of links related to this Commons Math GraggBulirschStoerStepInterpolator.java source code file:

Commons Math example source code file (GraggBulirschStoerStepInterpolator.java)

This example Commons Math source code file (GraggBulirschStoerStepInterpolator.java) is included in the DevDaily.com "Java Source Code Warehouse" project. The intent of this project is to help you "Learn Java by Example" TM.

Java - Commons Math tags/keywords

abstractstepinterpolator, derivativeexception, graggbulirschstoerstepinterpolator, graggbulirschstoerstepinterpolator, io, ioexception, ioexception, objectinput, objectinput, objectoutput, override, override, stepinterpolator

The Commons Math GraggBulirschStoerStepInterpolator.java source code

/*
 * 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.math.ode.nonstiff;

import java.io.IOException;
import java.io.ObjectInput;
import java.io.ObjectOutput;

import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
import org.apache.commons.math.ode.sampling.StepInterpolator;

/**
 * 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>
Copyright (c) 2004, Ernst Hairer
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>
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>
... this post is sponsored by my books ...

#1 New Release!

FP Best Seller

 

new blog posts

 

Copyright 1998-2024 Alvin Alexander, alvinalexander.com
All Rights Reserved.

A percentage of advertising revenue from
pages under the /java/jwarehouse URI on this website is
paid back to open source projects.