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