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