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