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