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 } 402
