Raymond | dee0849 | 2015-04-02 10:43:13 -0700 | [diff] [blame] | 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.sampling.AbstractStepInterpolator; |
| 25 | import org.apache.commons.math.ode.sampling.StepInterpolator; |
| 26 | import 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 | |
| 77 | class 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 | } |