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