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.math3.ode.nonstiff; 019 020 import org.apache.commons.math3.exception.DimensionMismatchException; 021 import org.apache.commons.math3.exception.MaxCountExceededException; 022 import org.apache.commons.math3.exception.NoBracketingException; 023 import org.apache.commons.math3.exception.NumberIsTooSmallException; 024 import org.apache.commons.math3.exception.util.LocalizedFormats; 025 import org.apache.commons.math3.ode.AbstractIntegrator; 026 import org.apache.commons.math3.ode.ExpandableStatefulODE; 027 import org.apache.commons.math3.util.FastMath; 028 029 /** 030 * This abstract class holds the common part of all adaptive 031 * stepsize integrators for Ordinary Differential Equations. 032 * 033 * <p>These algorithms perform integration with stepsize control, which 034 * means the user does not specify the integration step but rather a 035 * tolerance on error. The error threshold is computed as 036 * <pre> 037 * threshold_i = absTol_i + relTol_i * max (abs (ym), abs (ym+1)) 038 * </pre> 039 * where absTol_i is the absolute tolerance for component i of the 040 * state vector and relTol_i is the relative tolerance for the same 041 * component. The user can also use only two scalar values absTol and 042 * relTol which will be used for all components. 043 * </p> 044 * <p> 045 * If the Ordinary Differential Equations is an {@link ExpandableStatefulODE 046 * extended ODE} rather than a {@link 047 * org.apache.commons.math3.ode.FirstOrderDifferentialEquations basic ODE}, then 048 * <em>only</em> the {@link ExpandableStatefulODE#getPrimaryState() primary part} 049 * of the state vector is used for stepsize control, not the complete state vector. 050 * </p> 051 * 052 * <p>If the estimated error for ym+1 is such that 053 * <pre> 054 * sqrt((sum (errEst_i / threshold_i)^2 ) / n) < 1 055 * </pre> 056 * 057 * (where n is the main set dimension) then the step is accepted, 058 * otherwise the step is rejected and a new attempt is made with a new 059 * stepsize.</p> 060 * 061 * @version $Id: AdaptiveStepsizeIntegrator.java 1416643 2012-12-03 19:37:14Z tn $ 062 * @since 1.2 063 * 064 */ 065 066 public abstract class AdaptiveStepsizeIntegrator 067 extends AbstractIntegrator { 068 069 /** Allowed absolute scalar error. */ 070 protected double scalAbsoluteTolerance; 071 072 /** Allowed relative scalar error. */ 073 protected double scalRelativeTolerance; 074 075 /** Allowed absolute vectorial error. */ 076 protected double[] vecAbsoluteTolerance; 077 078 /** Allowed relative vectorial error. */ 079 protected double[] vecRelativeTolerance; 080 081 /** Main set dimension. */ 082 protected int mainSetDimension; 083 084 /** User supplied initial step. */ 085 private double initialStep; 086 087 /** Minimal step. */ 088 private double minStep; 089 090 /** Maximal step. */ 091 private double maxStep; 092 093 /** Build an integrator with the given stepsize bounds. 094 * The default step handler does nothing. 095 * @param name name of the method 096 * @param minStep minimal step (sign is irrelevant, regardless of 097 * integration direction, forward or backward), the last step can 098 * be smaller than this 099 * @param maxStep maximal step (sign is irrelevant, regardless of 100 * integration direction, forward or backward), the last step can 101 * be smaller than this 102 * @param scalAbsoluteTolerance allowed absolute error 103 * @param scalRelativeTolerance allowed relative error 104 */ 105 public AdaptiveStepsizeIntegrator(final String name, 106 final double minStep, final double maxStep, 107 final double scalAbsoluteTolerance, 108 final double scalRelativeTolerance) { 109 110 super(name); 111 setStepSizeControl(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); 112 resetInternalState(); 113 114 } 115 116 /** Build an integrator with the given stepsize bounds. 117 * The default step handler does nothing. 118 * @param name name of the method 119 * @param minStep minimal step (sign is irrelevant, regardless of 120 * integration direction, forward or backward), the last step can 121 * be smaller than this 122 * @param maxStep maximal step (sign is irrelevant, regardless of 123 * integration direction, forward or backward), the last step can 124 * be smaller than this 125 * @param vecAbsoluteTolerance allowed absolute error 126 * @param vecRelativeTolerance allowed relative error 127 */ 128 public AdaptiveStepsizeIntegrator(final String name, 129 final double minStep, final double maxStep, 130 final double[] vecAbsoluteTolerance, 131 final double[] vecRelativeTolerance) { 132 133 super(name); 134 setStepSizeControl(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); 135 resetInternalState(); 136 137 } 138 139 /** Set the adaptive step size control parameters. 140 * <p> 141 * A side effect of this method is to also reset the initial 142 * step so it will be automatically computed by the integrator 143 * if {@link #setInitialStepSize(double) setInitialStepSize} 144 * is not called by the user. 145 * </p> 146 * @param minimalStep minimal step (must be positive even for backward 147 * integration), the last step can be smaller than this 148 * @param maximalStep maximal step (must be positive even for backward 149 * integration) 150 * @param absoluteTolerance allowed absolute error 151 * @param relativeTolerance allowed relative error 152 */ 153 public void setStepSizeControl(final double minimalStep, final double maximalStep, 154 final double absoluteTolerance, 155 final double relativeTolerance) { 156 157 minStep = FastMath.abs(minimalStep); 158 maxStep = FastMath.abs(maximalStep); 159 initialStep = -1; 160 161 scalAbsoluteTolerance = absoluteTolerance; 162 scalRelativeTolerance = relativeTolerance; 163 vecAbsoluteTolerance = null; 164 vecRelativeTolerance = null; 165 166 } 167 168 /** Set the adaptive step size control parameters. 169 * <p> 170 * A side effect of this method is to also reset the initial 171 * step so it will be automatically computed by the integrator 172 * if {@link #setInitialStepSize(double) setInitialStepSize} 173 * is not called by the user. 174 * </p> 175 * @param minimalStep minimal step (must be positive even for backward 176 * integration), the last step can be smaller than this 177 * @param maximalStep maximal step (must be positive even for backward 178 * integration) 179 * @param absoluteTolerance allowed absolute error 180 * @param relativeTolerance allowed relative error 181 */ 182 public void setStepSizeControl(final double minimalStep, final double maximalStep, 183 final double[] absoluteTolerance, 184 final double[] relativeTolerance) { 185 186 minStep = FastMath.abs(minimalStep); 187 maxStep = FastMath.abs(maximalStep); 188 initialStep = -1; 189 190 scalAbsoluteTolerance = 0; 191 scalRelativeTolerance = 0; 192 vecAbsoluteTolerance = absoluteTolerance.clone(); 193 vecRelativeTolerance = relativeTolerance.clone(); 194 195 } 196 197 /** Set the initial step size. 198 * <p>This method allows the user to specify an initial positive 199 * step size instead of letting the integrator guess it by 200 * itself. If this method is not called before integration is 201 * started, the initial step size will be estimated by the 202 * integrator.</p> 203 * @param initialStepSize initial step size to use (must be positive even 204 * for backward integration ; providing a negative value or a value 205 * outside of the min/max step interval will lead the integrator to 206 * ignore the value and compute the initial step size by itself) 207 */ 208 public void setInitialStepSize(final double initialStepSize) { 209 if ((initialStepSize < minStep) || (initialStepSize > maxStep)) { 210 initialStep = -1.0; 211 } else { 212 initialStep = initialStepSize; 213 } 214 } 215 216 /** {@inheritDoc} */ 217 @Override 218 protected void sanityChecks(final ExpandableStatefulODE equations, final double t) 219 throws DimensionMismatchException, NumberIsTooSmallException { 220 221 super.sanityChecks(equations, t); 222 223 mainSetDimension = equations.getPrimaryMapper().getDimension(); 224 225 if ((vecAbsoluteTolerance != null) && (vecAbsoluteTolerance.length != mainSetDimension)) { 226 throw new DimensionMismatchException(mainSetDimension, vecAbsoluteTolerance.length); 227 } 228 229 if ((vecRelativeTolerance != null) && (vecRelativeTolerance.length != mainSetDimension)) { 230 throw new DimensionMismatchException(mainSetDimension, vecRelativeTolerance.length); 231 } 232 233 } 234 235 /** Initialize the integration step. 236 * @param forward forward integration indicator 237 * @param order order of the method 238 * @param scale scaling vector for the state vector (can be shorter than state vector) 239 * @param t0 start time 240 * @param y0 state vector at t0 241 * @param yDot0 first time derivative of y0 242 * @param y1 work array for a state vector 243 * @param yDot1 work array for the first time derivative of y1 244 * @return first integration step 245 * @exception MaxCountExceededException if the number of functions evaluations is exceeded 246 * @exception DimensionMismatchException if arrays dimensions do not match equations settings 247 */ 248 public double initializeStep(final boolean forward, final int order, final double[] scale, 249 final double t0, final double[] y0, final double[] yDot0, 250 final double[] y1, final double[] yDot1) 251 throws MaxCountExceededException, DimensionMismatchException { 252 253 if (initialStep > 0) { 254 // use the user provided value 255 return forward ? initialStep : -initialStep; 256 } 257 258 // very rough first guess : h = 0.01 * ||y/scale|| / ||y'/scale|| 259 // this guess will be used to perform an Euler step 260 double ratio; 261 double yOnScale2 = 0; 262 double yDotOnScale2 = 0; 263 for (int j = 0; j < scale.length; ++j) { 264 ratio = y0[j] / scale[j]; 265 yOnScale2 += ratio * ratio; 266 ratio = yDot0[j] / scale[j]; 267 yDotOnScale2 += ratio * ratio; 268 } 269 270 double h = ((yOnScale2 < 1.0e-10) || (yDotOnScale2 < 1.0e-10)) ? 271 1.0e-6 : (0.01 * FastMath.sqrt(yOnScale2 / yDotOnScale2)); 272 if (! forward) { 273 h = -h; 274 } 275 276 // perform an Euler step using the preceding rough guess 277 for (int j = 0; j < y0.length; ++j) { 278 y1[j] = y0[j] + h * yDot0[j]; 279 } 280 computeDerivatives(t0 + h, y1, yDot1); 281 282 // estimate the second derivative of the solution 283 double yDDotOnScale = 0; 284 for (int j = 0; j < scale.length; ++j) { 285 ratio = (yDot1[j] - yDot0[j]) / scale[j]; 286 yDDotOnScale += ratio * ratio; 287 } 288 yDDotOnScale = FastMath.sqrt(yDDotOnScale) / h; 289 290 // step size is computed such that 291 // h^order * max (||y'/tol||, ||y''/tol||) = 0.01 292 final double maxInv2 = FastMath.max(FastMath.sqrt(yDotOnScale2), yDDotOnScale); 293 final double h1 = (maxInv2 < 1.0e-15) ? 294 FastMath.max(1.0e-6, 0.001 * FastMath.abs(h)) : 295 FastMath.pow(0.01 / maxInv2, 1.0 / order); 296 h = FastMath.min(100.0 * FastMath.abs(h), h1); 297 h = FastMath.max(h, 1.0e-12 * FastMath.abs(t0)); // avoids cancellation when computing t1 - t0 298 if (h < getMinStep()) { 299 h = getMinStep(); 300 } 301 if (h > getMaxStep()) { 302 h = getMaxStep(); 303 } 304 if (! forward) { 305 h = -h; 306 } 307 308 return h; 309 310 } 311 312 /** Filter the integration step. 313 * @param h signed step 314 * @param forward forward integration indicator 315 * @param acceptSmall if true, steps smaller than the minimal value 316 * are silently increased up to this value, if false such small 317 * steps generate an exception 318 * @return a bounded integration step (h if no bound is reach, or a bounded value) 319 * @exception NumberIsTooSmallException if the step is too small and acceptSmall is false 320 */ 321 protected double filterStep(final double h, final boolean forward, final boolean acceptSmall) 322 throws NumberIsTooSmallException { 323 324 double filteredH = h; 325 if (FastMath.abs(h) < minStep) { 326 if (acceptSmall) { 327 filteredH = forward ? minStep : -minStep; 328 } else { 329 throw new NumberIsTooSmallException(LocalizedFormats.MINIMAL_STEPSIZE_REACHED_DURING_INTEGRATION, 330 FastMath.abs(h), minStep, true); 331 } 332 } 333 334 if (filteredH > maxStep) { 335 filteredH = maxStep; 336 } else if (filteredH < -maxStep) { 337 filteredH = -maxStep; 338 } 339 340 return filteredH; 341 342 } 343 344 /** {@inheritDoc} */ 345 @Override 346 public abstract void integrate (ExpandableStatefulODE equations, double t) 347 throws NumberIsTooSmallException, DimensionMismatchException, 348 MaxCountExceededException, NoBracketingException; 349 350 /** {@inheritDoc} */ 351 @Override 352 public double getCurrentStepStart() { 353 return stepStart; 354 } 355 356 /** Reset internal state to dummy values. */ 357 protected void resetInternalState() { 358 stepStart = Double.NaN; 359 stepSize = FastMath.sqrt(minStep * maxStep); 360 } 361 362 /** Get the minimal step. 363 * @return minimal step 364 */ 365 public double getMinStep() { 366 return minStep; 367 } 368 369 /** Get the maximal step. 370 * @return maximal step 371 */ 372 public double getMaxStep() { 373 return maxStep; 374 } 375 376 }