001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018 package org.apache.commons.math.ode; 019 020 import org.apache.commons.math.MathRuntimeException; 021 import org.apache.commons.math.linear.Array2DRowRealMatrix; 022 import org.apache.commons.math.linear.RealMatrix; 023 import org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator; 024 import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator; 025 import org.apache.commons.math.ode.sampling.StepHandler; 026 import org.apache.commons.math.ode.sampling.StepInterpolator; 027 028 /** 029 * This class is the base class for multistep integrators for Ordinary 030 * Differential Equations. 031 * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as: 032 * <pre> 033 * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative 034 * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative 035 * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative 036 * ... 037 * s<sub>k</sub>(n) = h<sup>k</sup>/k! y(k)<sub>n</sub> for k<sup>th</sup> derivative 038 * </pre></p> 039 * <p>Rather than storing several previous steps separately, this implementation uses 040 * the Nordsieck vector with higher degrees scaled derivatives all taken at the same 041 * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as: 042 * <pre> 043 * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup> 044 * </pre> 045 * (we omit the k index in the notation for clarity)</p> 046 * <p> 047 * Multistep integrators with Nordsieck representation are highly sensitive to 048 * large step changes because when the step is multiplied by a factor a, the 049 * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup> 050 * and the last components are the least accurate ones. The default max growth 051 * factor is therefore set to a quite low value: 2<sup>1/order</sup>. 052 * </p> 053 * 054 * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator 055 * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator 056 * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $ 057 * @since 2.0 058 */ 059 public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator { 060 061 /** Starter integrator. */ 062 private FirstOrderIntegrator starter; 063 064 /** Number of steps of the multistep method (excluding the one being computed). */ 065 private final int nSteps; 066 067 /** First scaled derivative (h y'). */ 068 protected double[] scaled; 069 070 /** Nordsieck matrix of the higher scaled derivatives. 071 * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y(k))</p> 072 */ 073 protected Array2DRowRealMatrix nordsieck; 074 075 /** Stepsize control exponent. */ 076 private double exp; 077 078 /** Safety factor for stepsize control. */ 079 private double safety; 080 081 /** Minimal reduction factor for stepsize control. */ 082 private double minReduction; 083 084 /** Maximal growth factor for stepsize control. */ 085 private double maxGrowth; 086 087 /** 088 * Build a multistep integrator with the given stepsize bounds. 089 * <p>The default starter integrator is set to the {@link 090 * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with 091 * some defaults settings.</p> 092 * <p> 093 * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>. 094 * </p> 095 * @param name name of the method 096 * @param nSteps number of steps of the multistep method 097 * (excluding the one being computed) 098 * @param order order of the method 099 * @param minStep minimal step (must be positive even for backward 100 * integration), the last step can be smaller than this 101 * @param maxStep maximal step (must be positive even for backward 102 * integration) 103 * @param scalAbsoluteTolerance allowed absolute error 104 * @param scalRelativeTolerance allowed relative error 105 */ 106 protected MultistepIntegrator(final String name, final int nSteps, 107 final int order, 108 final double minStep, final double maxStep, 109 final double scalAbsoluteTolerance, 110 final double scalRelativeTolerance) { 111 112 super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); 113 114 if (nSteps <= 0) { 115 throw MathRuntimeException.createIllegalArgumentException( 116 "{0} method needs at least one previous point", 117 name); 118 } 119 120 starter = new DormandPrince853Integrator(minStep, maxStep, 121 scalAbsoluteTolerance, 122 scalRelativeTolerance); 123 this.nSteps = nSteps; 124 125 exp = -1.0 / order; 126 127 // set the default values of the algorithm control parameters 128 setSafety(0.9); 129 setMinReduction(0.2); 130 setMaxGrowth(Math.pow(2.0, -exp)); 131 132 } 133 134 /** 135 * Build a multistep integrator with the given stepsize bounds. 136 * <p>The default starter integrator is set to the {@link 137 * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with 138 * some defaults settings.</p> 139 * <p> 140 * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>. 141 * </p> 142 * @param name name of the method 143 * @param nSteps number of steps of the multistep method 144 * (excluding the one being computed) 145 * @param order order of the method 146 * @param minStep minimal step (must be positive even for backward 147 * integration), the last step can be smaller than this 148 * @param maxStep maximal step (must be positive even for backward 149 * integration) 150 * @param vecAbsoluteTolerance allowed absolute error 151 * @param vecRelativeTolerance allowed relative error 152 */ 153 protected MultistepIntegrator(final String name, final int nSteps, 154 final int order, 155 final double minStep, final double maxStep, 156 final double[] vecAbsoluteTolerance, 157 final double[] vecRelativeTolerance) { 158 super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); 159 starter = new DormandPrince853Integrator(minStep, maxStep, 160 vecAbsoluteTolerance, 161 vecRelativeTolerance); 162 this.nSteps = nSteps; 163 164 exp = -1.0 / order; 165 166 // set the default values of the algorithm control parameters 167 setSafety(0.9); 168 setMinReduction(0.2); 169 setMaxGrowth(Math.pow(2.0, -exp)); 170 171 } 172 173 /** 174 * Get the starter integrator. 175 * @return starter integrator 176 */ 177 public ODEIntegrator getStarterIntegrator() { 178 return starter; 179 } 180 181 /** 182 * Set the starter integrator. 183 * <p>The various step and event handlers for this starter integrator 184 * will be managed automatically by the multi-step integrator. Any 185 * user configuration for these elements will be cleared before use.</p> 186 * @param starter starter integrator 187 */ 188 public void setStarterIntegrator(FirstOrderIntegrator starter) { 189 this.starter = starter; 190 } 191 192 /** Start the integration. 193 * <p>This method computes one step using the underlying starter integrator, 194 * and initializes the Nordsieck vector at step start. The starter integrator 195 * purpose is only to establish initial conditions, it does not really change 196 * time by itself. The top level multistep integrator remains in charge of 197 * handling time propagation and events handling as it will starts its own 198 * computation right from the beginning. In a sense, the starter integrator 199 * can be seen as a dummy one and so it will never trigger any user event nor 200 * call any user step handler.</p> 201 * @param t0 initial time 202 * @param y0 initial value of the state vector at t0 203 * @param t target time for the integration 204 * (can be set to a value smaller than <code>t0</code> for backward integration) 205 * @throws IntegratorException if the integrator cannot perform integration 206 * @throws DerivativeException this exception is propagated to the caller if 207 * the underlying user function triggers one 208 */ 209 protected void start(final double t0, final double[] y0, final double t) 210 throws DerivativeException, IntegratorException { 211 212 // make sure NO user event nor user step handler is triggered, 213 // this is the task of the top level integrator, not the task 214 // of the starter integrator 215 starter.clearEventHandlers(); 216 starter.clearStepHandlers(); 217 218 // set up one specific step handler to extract initial Nordsieck vector 219 starter.addStepHandler(new NordsieckInitializer(y0.length)); 220 221 // start integration, expecting a InitializationCompletedMarkerException 222 try { 223 starter.integrate(new CountingDifferentialEquations(y0.length), 224 t0, y0, t, new double[y0.length]); 225 } catch (DerivativeException de) { 226 if (!(de instanceof InitializationCompletedMarkerException)) { 227 // this is not the expected nominal interruption of the start integrator 228 throw de; 229 } 230 } 231 232 // remove the specific step handler 233 starter.clearStepHandlers(); 234 235 } 236 237 /** Initialize the high order scaled derivatives at step start. 238 * @param first first scaled derivative at step start 239 * @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1) 240 * will be modified 241 * @return high order scaled derivatives at step start 242 */ 243 protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first, 244 final double[][] multistep); 245 246 /** Get the minimal reduction factor for stepsize control. 247 * @return minimal reduction factor 248 */ 249 public double getMinReduction() { 250 return minReduction; 251 } 252 253 /** Set the minimal reduction factor for stepsize control. 254 * @param minReduction minimal reduction factor 255 */ 256 public void setMinReduction(final double minReduction) { 257 this.minReduction = minReduction; 258 } 259 260 /** Get the maximal growth factor for stepsize control. 261 * @return maximal growth factor 262 */ 263 public double getMaxGrowth() { 264 return maxGrowth; 265 } 266 267 /** Set the maximal growth factor for stepsize control. 268 * @param maxGrowth maximal growth factor 269 */ 270 public void setMaxGrowth(final double maxGrowth) { 271 this.maxGrowth = maxGrowth; 272 } 273 274 /** Get the safety factor for stepsize control. 275 * @return safety factor 276 */ 277 public double getSafety() { 278 return safety; 279 } 280 281 /** Set the safety factor for stepsize control. 282 * @param safety safety factor 283 */ 284 public void setSafety(final double safety) { 285 this.safety = safety; 286 } 287 288 /** Compute step grow/shrink factor according to normalized error. 289 * @param error normalized error of the current step 290 * @return grow/shrink factor for next step 291 */ 292 protected double computeStepGrowShrinkFactor(final double error) { 293 return Math.min(maxGrowth, Math.max(minReduction, safety * Math.pow(error, exp))); 294 } 295 296 /** Transformer used to convert the first step to Nordsieck representation. */ 297 public static interface NordsieckTransformer { 298 /** Initialize the high order scaled derivatives at step start. 299 * @param first first scaled derivative at step start 300 * @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1) 301 * will be modified 302 * @return high order derivatives at step start 303 */ 304 RealMatrix initializeHighOrderDerivatives(double[] first, double[][] multistep); 305 } 306 307 /** Specialized step handler storing the first step. */ 308 private class NordsieckInitializer implements StepHandler { 309 310 /** Problem dimension. */ 311 private final int n; 312 313 /** Simple constructor. 314 * @param n problem dimension 315 */ 316 public NordsieckInitializer(final int n) { 317 this.n = n; 318 } 319 320 /** {@inheritDoc} */ 321 public void handleStep(StepInterpolator interpolator, boolean isLast) 322 throws DerivativeException { 323 324 final double prev = interpolator.getPreviousTime(); 325 final double curr = interpolator.getCurrentTime(); 326 stepStart = prev; 327 stepSize = (curr - prev) / (nSteps + 1); 328 329 // compute the first scaled derivative 330 interpolator.setInterpolatedTime(prev); 331 scaled = interpolator.getInterpolatedDerivatives().clone(); 332 for (int j = 0; j < n; ++j) { 333 scaled[j] *= stepSize; 334 } 335 336 // compute the high order scaled derivatives 337 final double[][] multistep = new double[nSteps][]; 338 for (int i = 1; i <= nSteps; ++i) { 339 interpolator.setInterpolatedTime(prev + stepSize * i); 340 final double[] msI = interpolator.getInterpolatedDerivatives().clone(); 341 for (int j = 0; j < n; ++j) { 342 msI[j] *= stepSize; 343 } 344 multistep[i - 1] = msI; 345 } 346 nordsieck = initializeHighOrderDerivatives(scaled, multistep); 347 348 // stop the integrator after the first step has been handled 349 throw new InitializationCompletedMarkerException(); 350 351 } 352 353 /** {@inheritDoc} */ 354 public boolean requiresDenseOutput() { 355 return true; 356 } 357 358 /** {@inheritDoc} */ 359 public void reset() { 360 // nothing to do 361 } 362 363 } 364 365 /** Marker exception used ONLY to stop the starter integrator after first step. */ 366 private static class InitializationCompletedMarkerException 367 extends DerivativeException { 368 369 /** Serializable version identifier. */ 370 private static final long serialVersionUID = -4105805787353488365L; 371 372 /** Simple constructor. */ 373 public InitializationCompletedMarkerException() { 374 super((Throwable) null); 375 } 376 377 } 378 379 /** Wrapper for differential equations, ensuring start evaluations are counted. */ 380 private class CountingDifferentialEquations implements FirstOrderDifferentialEquations { 381 382 /** Dimension of the problem. */ 383 private final int dimension; 384 385 /** Simple constructor. 386 * @param dimension dimension of the problem 387 */ 388 public CountingDifferentialEquations(final int dimension) { 389 this.dimension = dimension; 390 } 391 392 /** {@inheritDoc} */ 393 public void computeDerivatives(double t, double[] y, double[] dot) 394 throws DerivativeException { 395 MultistepIntegrator.this.computeDerivatives(t, y, dot); 396 } 397 398 /** {@inheritDoc} */ 399 public int getDimension() { 400 return dimension; 401 } 402 } 403 404 }