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