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