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