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