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