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; 19 20import org.apache.commons.math.MathRuntimeException; 21import org.apache.commons.math.ode.DerivativeException; 22import org.apache.commons.math.exception.util.LocalizedFormats; 23import org.apache.commons.math.linear.Array2DRowRealMatrix; 24import org.apache.commons.math.linear.RealMatrix; 25import org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator; 26import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator; 27import org.apache.commons.math.ode.sampling.StepHandler; 28import org.apache.commons.math.ode.sampling.StepInterpolator; 29import org.apache.commons.math.util.FastMath; 30 31/** 32 * This class is the base class for multistep integrators for Ordinary 33 * Differential Equations. 34 * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as: 35 * <pre> 36 * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative 37 * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative 38 * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative 39 * ... 40 * s<sub>k</sub>(n) = h<sup>k</sup>/k! y(k)<sub>n</sub> for k<sup>th</sup> derivative 41 * </pre></p> 42 * <p>Rather than storing several previous steps separately, this implementation uses 43 * the Nordsieck vector with higher degrees scaled derivatives all taken at the same 44 * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as: 45 * <pre> 46 * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup> 47 * </pre> 48 * (we omit the k index in the notation for clarity)</p> 49 * <p> 50 * Multistep integrators with Nordsieck representation are highly sensitive to 51 * large step changes because when the step is multiplied by a factor a, the 52 * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup> 53 * and the last components are the least accurate ones. The default max growth 54 * factor is therefore set to a quite low value: 2<sup>1/order</sup>. 55 * </p> 56 * 57 * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator 58 * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator 59 * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $ 60 * @since 2.0 61 */ 62public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator { 63 64 /** First scaled derivative (h y'). */ 65 protected double[] scaled; 66 67 /** Nordsieck matrix of the higher scaled derivatives. 68 * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y(k))</p> 69 */ 70 protected Array2DRowRealMatrix nordsieck; 71 72 /** Starter integrator. */ 73 private FirstOrderIntegrator starter; 74 75 /** Number of steps of the multistep method (excluding the one being computed). */ 76 private final int nSteps; 77 78 /** Stepsize control exponent. */ 79 private double exp; 80 81 /** Safety factor for stepsize control. */ 82 private double safety; 83 84 /** Minimal reduction factor for stepsize control. */ 85 private double minReduction; 86 87 /** Maximal growth factor for stepsize control. */ 88 private double maxGrowth; 89 90 /** 91 * Build a multistep integrator with the given stepsize bounds. 92 * <p>The default starter integrator is set to the {@link 93 * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with 94 * some defaults settings.</p> 95 * <p> 96 * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>. 97 * </p> 98 * @param name name of the method 99 * @param nSteps number of steps of the multistep method 100 * (excluding the one being computed) 101 * @param order order of the method 102 * @param minStep minimal step (must be positive even for backward 103 * integration), the last step can be smaller than this 104 * @param maxStep maximal step (must be positive even for backward 105 * integration) 106 * @param scalAbsoluteTolerance allowed absolute error 107 * @param scalRelativeTolerance allowed relative error 108 */ 109 protected MultistepIntegrator(final String name, final int nSteps, 110 final int order, 111 final double minStep, final double maxStep, 112 final double scalAbsoluteTolerance, 113 final double scalRelativeTolerance) { 114 115 super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); 116 117 if (nSteps <= 0) { 118 throw MathRuntimeException.createIllegalArgumentException( 119 LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_ONE_PREVIOUS_POINT, 120 name); 121 } 122 123 starter = new DormandPrince853Integrator(minStep, maxStep, 124 scalAbsoluteTolerance, 125 scalRelativeTolerance); 126 this.nSteps = nSteps; 127 128 exp = -1.0 / order; 129 130 // set the default values of the algorithm control parameters 131 setSafety(0.9); 132 setMinReduction(0.2); 133 setMaxGrowth(FastMath.pow(2.0, -exp)); 134 135 } 136 137 /** 138 * Build a multistep integrator with the given stepsize bounds. 139 * <p>The default starter integrator is set to the {@link 140 * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with 141 * some defaults settings.</p> 142 * <p> 143 * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>. 144 * </p> 145 * @param name name of the method 146 * @param nSteps number of steps of the multistep method 147 * (excluding the one being computed) 148 * @param order order of the method 149 * @param minStep minimal step (must be positive even for backward 150 * integration), the last step can be smaller than this 151 * @param maxStep maximal step (must be positive even for backward 152 * integration) 153 * @param vecAbsoluteTolerance allowed absolute error 154 * @param vecRelativeTolerance allowed relative error 155 */ 156 protected MultistepIntegrator(final String name, final int nSteps, 157 final int order, 158 final double minStep, final double maxStep, 159 final double[] vecAbsoluteTolerance, 160 final double[] vecRelativeTolerance) { 161 super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); 162 starter = new DormandPrince853Integrator(minStep, maxStep, 163 vecAbsoluteTolerance, 164 vecRelativeTolerance); 165 this.nSteps = nSteps; 166 167 exp = -1.0 / order; 168 169 // set the default values of the algorithm control parameters 170 setSafety(0.9); 171 setMinReduction(0.2); 172 setMaxGrowth(FastMath.pow(2.0, -exp)); 173 174 } 175 176 /** 177 * Get the starter integrator. 178 * @return starter integrator 179 */ 180 public ODEIntegrator getStarterIntegrator() { 181 return starter; 182 } 183 184 /** 185 * Set the starter integrator. 186 * <p>The various step and event handlers for this starter integrator 187 * will be managed automatically by the multi-step integrator. Any 188 * user configuration for these elements will be cleared before use.</p> 189 * @param starterIntegrator starter integrator 190 */ 191 public void setStarterIntegrator(FirstOrderIntegrator starterIntegrator) { 192 this.starter = starterIntegrator; 193 } 194 195 /** Start the integration. 196 * <p>This method computes one step using the underlying starter integrator, 197 * and initializes the Nordsieck vector at step start. The starter integrator 198 * purpose is only to establish initial conditions, it does not really change 199 * time by itself. The top level multistep integrator remains in charge of 200 * handling time propagation and events handling as it will starts its own 201 * computation right from the beginning. In a sense, the starter integrator 202 * can be seen as a dummy one and so it will never trigger any user event nor 203 * call any user step handler.</p> 204 * @param t0 initial time 205 * @param y0 initial value of the state vector at t0 206 * @param t target time for the integration 207 * (can be set to a value smaller than <code>t0</code> for backward integration) 208 * @throws IntegratorException if the integrator cannot perform integration 209 * @throws DerivativeException this exception is propagated to the caller if 210 * the underlying user function triggers one 211 */ 212 protected void start(final double t0, final double[] y0, final double t) 213 throws DerivativeException, IntegratorException { 214 215 // make sure NO user event nor user step handler is triggered, 216 // this is the task of the top level integrator, not the task 217 // of the starter integrator 218 starter.clearEventHandlers(); 219 starter.clearStepHandlers(); 220 221 // set up one specific step handler to extract initial Nordsieck vector 222 starter.addStepHandler(new NordsieckInitializer(y0.length)); 223 224 // start integration, expecting a InitializationCompletedMarkerException 225 try { 226 starter.integrate(new CountingDifferentialEquations(y0.length), 227 t0, y0, t, new double[y0.length]); 228 } catch (DerivativeException mue) { 229 if (!(mue instanceof InitializationCompletedMarkerException)) { 230 // this is not the expected nominal interruption of the start integrator 231 throw mue; 232 } 233 } 234 235 // remove the specific step handler 236 starter.clearStepHandlers(); 237 238 } 239 240 /** Initialize the high order scaled derivatives at step start. 241 * @param first first scaled derivative at step start 242 * @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1) 243 * will be modified 244 * @return high order scaled derivatives at step start 245 */ 246 protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first, 247 final double[][] multistep); 248 249 /** Get the minimal reduction factor for stepsize control. 250 * @return minimal reduction factor 251 */ 252 public double getMinReduction() { 253 return minReduction; 254 } 255 256 /** Set the minimal reduction factor for stepsize control. 257 * @param minReduction minimal reduction factor 258 */ 259 public void setMinReduction(final double minReduction) { 260 this.minReduction = minReduction; 261 } 262 263 /** Get the maximal growth factor for stepsize control. 264 * @return maximal growth factor 265 */ 266 public double getMaxGrowth() { 267 return maxGrowth; 268 } 269 270 /** Set the maximal growth factor for stepsize control. 271 * @param maxGrowth maximal growth factor 272 */ 273 public void setMaxGrowth(final double maxGrowth) { 274 this.maxGrowth = maxGrowth; 275 } 276 277 /** Get the safety factor for stepsize control. 278 * @return safety factor 279 */ 280 public double getSafety() { 281 return safety; 282 } 283 284 /** Set the safety factor for stepsize control. 285 * @param safety safety factor 286 */ 287 public void setSafety(final double safety) { 288 this.safety = safety; 289 } 290 291 /** Compute step grow/shrink factor according to normalized error. 292 * @param error normalized error of the current step 293 * @return grow/shrink factor for next step 294 */ 295 protected double computeStepGrowShrinkFactor(final double error) { 296 return FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp))); 297 } 298 299 /** Transformer used to convert the first step to Nordsieck representation. */ 300 public static interface NordsieckTransformer { 301 /** Initialize the high order scaled derivatives at step start. 302 * @param first first scaled derivative at step start 303 * @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1) 304 * will be modified 305 * @return high order derivatives at step start 306 */ 307 RealMatrix initializeHighOrderDerivatives(double[] first, double[][] multistep); 308 } 309 310 /** Specialized step handler storing the first step. */ 311 private class NordsieckInitializer implements StepHandler { 312 313 /** Problem dimension. */ 314 private final int n; 315 316 /** Simple constructor. 317 * @param n problem dimension 318 */ 319 public NordsieckInitializer(final int n) { 320 this.n = n; 321 } 322 323 /** {@inheritDoc} */ 324 public void handleStep(StepInterpolator interpolator, boolean isLast) 325 throws DerivativeException { 326 327 final double prev = interpolator.getPreviousTime(); 328 final double curr = interpolator.getCurrentTime(); 329 stepStart = prev; 330 stepSize = (curr - prev) / (nSteps + 1); 331 332 // compute the first scaled derivative 333 interpolator.setInterpolatedTime(prev); 334 scaled = interpolator.getInterpolatedDerivatives().clone(); 335 for (int j = 0; j < n; ++j) { 336 scaled[j] *= stepSize; 337 } 338 339 // compute the high order scaled derivatives 340 final double[][] multistep = new double[nSteps][]; 341 for (int i = 1; i <= nSteps; ++i) { 342 interpolator.setInterpolatedTime(prev + stepSize * i); 343 final double[] msI = interpolator.getInterpolatedDerivatives().clone(); 344 for (int j = 0; j < n; ++j) { 345 msI[j] *= stepSize; 346 } 347 multistep[i - 1] = msI; 348 } 349 nordsieck = initializeHighOrderDerivatives(scaled, multistep); 350 351 // stop the integrator after the first step has been handled 352 throw new InitializationCompletedMarkerException(); 353 354 } 355 356 /** {@inheritDoc} */ 357 public boolean requiresDenseOutput() { 358 return true; 359 } 360 361 /** {@inheritDoc} */ 362 public void reset() { 363 // nothing to do 364 } 365 366 } 367 368 /** Marker exception used ONLY to stop the starter integrator after first step. */ 369 private static class InitializationCompletedMarkerException 370 extends DerivativeException { 371 372 /** Serializable version identifier. */ 373 private static final long serialVersionUID = -4105805787353488365L; 374 375 /** Simple constructor. */ 376 public InitializationCompletedMarkerException() { 377 super((Throwable) null); 378 } 379 380 } 381 382 /** Wrapper for differential equations, ensuring start evaluations are counted. */ 383 private class CountingDifferentialEquations implements ExtendedFirstOrderDifferentialEquations { 384 385 /** Dimension of the problem. */ 386 private final int dimension; 387 388 /** Simple constructor. 389 * @param dimension dimension of the problem 390 */ 391 public CountingDifferentialEquations(final int dimension) { 392 this.dimension = dimension; 393 } 394 395 /** {@inheritDoc} */ 396 public void computeDerivatives(double t, double[] y, double[] dot) 397 throws DerivativeException { 398 MultistepIntegrator.this.computeDerivatives(t, y, dot); 399 } 400 401 /** {@inheritDoc} */ 402 public int getDimension() { 403 return dimension; 404 } 405 406 /** {@inheritDoc} */ 407 public int getMainSetDimension() { 408 return mainSetDimension; 409 } 410 } 411 412} 413