1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 package org.apache.commons.math4.legacy.ode.nonstiff;
19
20 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
21 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
22 import org.apache.commons.math4.legacy.exception.NoBracketingException;
23 import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
24 import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
25 import org.apache.commons.math4.legacy.ode.AbstractIntegrator;
26 import org.apache.commons.math4.legacy.ode.ExpandableStatefulODE;
27 import org.apache.commons.math4.core.jdkmath.JdkMath;
28
29 /**
30 * This abstract class holds the common part of all adaptive
31 * stepsize integrators for Ordinary Differential Equations.
32 *
33 * <p>These algorithms perform integration with stepsize control, which
34 * means the user does not specify the integration step but rather a
35 * tolerance on error. The error threshold is computed as
36 * <pre>
37 * threshold_i = absTol_i + relTol_i * max (abs (ym), abs (ym+1))
38 * </pre>
39 * where absTol_i is the absolute tolerance for component i of the
40 * state vector and relTol_i is the relative tolerance for the same
41 * component. The user can also use only two scalar values absTol and
42 * relTol which will be used for all components.
43 *
44 * <p>
45 * If the Ordinary Differential Equations is an {@link ExpandableStatefulODE
46 * extended ODE} rather than a {@link
47 * org.apache.commons.math4.legacy.ode.FirstOrderDifferentialEquations basic ODE}, then
48 * <em>only</em> the {@link ExpandableStatefulODE#getPrimaryState() primary part}
49 * of the state vector is used for stepsize control, not the complete state vector.
50 * </p>
51 *
52 * <p>If the estimated error for ym+1 is such that
53 * <pre>
54 * sqrt((sum (errEst_i / threshold_i)^2 ) / n) < 1
55 * </pre>
56 *
57 * (where n is the main set dimension) then the step is accepted,
58 * otherwise the step is rejected and a new attempt is made with a new
59 * stepsize.
60 *
61 * @since 1.2
62 *
63 */
64
65 public abstract class AdaptiveStepsizeIntegrator
66 extends AbstractIntegrator {
67
68 /** Allowed absolute scalar error. */
69 protected double scalAbsoluteTolerance;
70
71 /** Allowed relative scalar error. */
72 protected double scalRelativeTolerance;
73
74 /** Allowed absolute vectorial error. */
75 protected double[] vecAbsoluteTolerance;
76
77 /** Allowed relative vectorial error. */
78 protected double[] vecRelativeTolerance;
79
80 /** Main set dimension. */
81 protected int mainSetDimension;
82
83 /** User supplied initial step. */
84 private double initialStep;
85
86 /** Minimal step. */
87 private double minStep;
88
89 /** Maximal step. */
90 private double maxStep;
91
92 /** Build an integrator with the given stepsize bounds.
93 * The default step handler does nothing.
94 * @param name name of the method
95 * @param minStep minimal step (sign is irrelevant, regardless of
96 * integration direction, forward or backward), the last step can
97 * be smaller than this
98 * @param maxStep maximal step (sign is irrelevant, regardless of
99 * integration direction, forward or backward), the last step can
100 * be smaller than this
101 * @param scalAbsoluteTolerance allowed absolute error
102 * @param scalRelativeTolerance allowed relative error
103 */
104 public AdaptiveStepsizeIntegrator(final String name,
105 final double minStep, final double maxStep,
106 final double scalAbsoluteTolerance,
107 final double scalRelativeTolerance) {
108
109 super(name);
110 setStepSizeControl(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
111 resetInternalState();
112 }
113
114 /** Build an integrator with the given stepsize bounds.
115 * The default step handler does nothing.
116 * @param name name of the method
117 * @param minStep minimal step (sign is irrelevant, regardless of
118 * integration direction, forward or backward), the last step can
119 * be smaller than this
120 * @param maxStep maximal step (sign is irrelevant, regardless of
121 * integration direction, forward or backward), the last step can
122 * be smaller than this
123 * @param vecAbsoluteTolerance allowed absolute error
124 * @param vecRelativeTolerance allowed relative error
125 */
126 public AdaptiveStepsizeIntegrator(final String name,
127 final double minStep, final double maxStep,
128 final double[] vecAbsoluteTolerance,
129 final double[] vecRelativeTolerance) {
130
131 super(name);
132 setStepSizeControl(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
133 resetInternalState();
134 }
135
136 /** Set the adaptive step size control parameters.
137 * <p>
138 * A side effect of this method is to also reset the initial
139 * step so it will be automatically computed by the integrator
140 * if {@link #setInitialStepSize(double) setInitialStepSize}
141 * is not called by the user.
142 * </p>
143 * @param minimalStep minimal step (must be positive even for backward
144 * integration), the last step can be smaller than this
145 * @param maximalStep maximal step (must be positive even for backward
146 * integration)
147 * @param absoluteTolerance allowed absolute error
148 * @param relativeTolerance allowed relative error
149 */
150 public void setStepSizeControl(final double minimalStep, final double maximalStep,
151 final double absoluteTolerance,
152 final double relativeTolerance) {
153
154 minStep = JdkMath.abs(minimalStep);
155 maxStep = JdkMath.abs(maximalStep);
156 initialStep = -1;
157
158 scalAbsoluteTolerance = absoluteTolerance;
159 scalRelativeTolerance = relativeTolerance;
160 vecAbsoluteTolerance = null;
161 vecRelativeTolerance = null;
162 }
163
164 /** Set the adaptive step size control parameters.
165 * <p>
166 * A side effect of this method is to also reset the initial
167 * step so it will be automatically computed by the integrator
168 * if {@link #setInitialStepSize(double) setInitialStepSize}
169 * is not called by the user.
170 * </p>
171 * @param minimalStep minimal step (must be positive even for backward
172 * integration), the last step can be smaller than this
173 * @param maximalStep maximal step (must be positive even for backward
174 * integration)
175 * @param absoluteTolerance allowed absolute error
176 * @param relativeTolerance allowed relative error
177 */
178 public void setStepSizeControl(final double minimalStep, final double maximalStep,
179 final double[] absoluteTolerance,
180 final double[] relativeTolerance) {
181
182 minStep = JdkMath.abs(minimalStep);
183 maxStep = JdkMath.abs(maximalStep);
184 initialStep = -1;
185
186 scalAbsoluteTolerance = 0;
187 scalRelativeTolerance = 0;
188 vecAbsoluteTolerance = absoluteTolerance.clone();
189 vecRelativeTolerance = relativeTolerance.clone();
190 }
191
192 /** Set the initial step size.
193 * <p>This method allows the user to specify an initial positive
194 * step size instead of letting the integrator guess it by
195 * itself. If this method is not called before integration is
196 * started, the initial step size will be estimated by the
197 * integrator.</p>
198 * @param initialStepSize initial step size to use (must be positive even
199 * for backward integration ; providing a negative value or a value
200 * outside of the min/max step interval will lead the integrator to
201 * ignore the value and compute the initial step size by itself)
202 */
203 public void setInitialStepSize(final double initialStepSize) {
204 if (initialStepSize < minStep || initialStepSize > maxStep) {
205 initialStep = -1.0;
206 } else {
207 initialStep = initialStepSize;
208 }
209 }
210
211 /** {@inheritDoc} */
212 @Override
213 protected void sanityChecks(final ExpandableStatefulODE equations, final double t)
214 throws DimensionMismatchException, NumberIsTooSmallException {
215
216 super.sanityChecks(equations, t);
217
218 mainSetDimension = equations.getPrimaryMapper().getDimension();
219
220 if (vecAbsoluteTolerance != null && vecAbsoluteTolerance.length != mainSetDimension) {
221 throw new DimensionMismatchException(mainSetDimension, vecAbsoluteTolerance.length);
222 }
223
224 if (vecRelativeTolerance != null && vecRelativeTolerance.length != mainSetDimension) {
225 throw new DimensionMismatchException(mainSetDimension, vecRelativeTolerance.length);
226 }
227 }
228
229 /** Initialize the integration step.
230 * @param forward forward integration indicator
231 * @param order order of the method
232 * @param scale scaling vector for the state vector (can be shorter than state vector)
233 * @param t0 start time
234 * @param y0 state vector at t0
235 * @param yDot0 first time derivative of y0
236 * @param y1 work array for a state vector
237 * @param yDot1 work array for the first time derivative of y1
238 * @return first integration step
239 * @exception MaxCountExceededException if the number of functions evaluations is exceeded
240 * @exception DimensionMismatchException if arrays dimensions do not match equations settings
241 */
242 public double initializeStep(final boolean forward, final int order, final double[] scale,
243 final double t0, final double[] y0, final double[] yDot0,
244 final double[] y1, final double[] yDot1)
245 throws MaxCountExceededException, DimensionMismatchException {
246
247 if (initialStep > 0) {
248 // use the user provided value
249 return forward ? initialStep : -initialStep;
250 }
251
252 // very rough first guess : h = 0.01 * ||y/scale|| / ||y'/scale||
253 // this guess will be used to perform an Euler step
254 double ratio;
255 double yOnScale2 = 0;
256 double yDotOnScale2 = 0;
257 for (int j = 0; j < scale.length; ++j) {
258 ratio = y0[j] / scale[j];
259 yOnScale2 += ratio * ratio;
260 ratio = yDot0[j] / scale[j];
261 yDotOnScale2 += ratio * ratio;
262 }
263
264 double h = (yOnScale2 < 1.0e-10 || yDotOnScale2 < 1.0e-10) ?
265 1.0e-6 : (0.01 * JdkMath.sqrt(yOnScale2 / yDotOnScale2));
266 if (! forward) {
267 h = -h;
268 }
269
270 // perform an Euler step using the preceding rough guess
271 for (int j = 0; j < y0.length; ++j) {
272 y1[j] = y0[j] + h * yDot0[j];
273 }
274 computeDerivatives(t0 + h, y1, yDot1);
275
276 // estimate the second derivative of the solution
277 double yDDotOnScale = 0;
278 for (int j = 0; j < scale.length; ++j) {
279 ratio = (yDot1[j] - yDot0[j]) / scale[j];
280 yDDotOnScale += ratio * ratio;
281 }
282 yDDotOnScale = JdkMath.sqrt(yDDotOnScale) / h;
283
284 // step size is computed such that
285 // h^order * max (||y'/tol||, ||y''/tol||) = 0.01
286 final double maxInv2 = JdkMath.max(JdkMath.sqrt(yDotOnScale2), yDDotOnScale);
287 final double h1 = (maxInv2 < 1.0e-15) ?
288 JdkMath.max(1.0e-6, 0.001 * JdkMath.abs(h)) :
289 JdkMath.pow(0.01 / maxInv2, 1.0 / order);
290 h = JdkMath.min(100.0 * JdkMath.abs(h), h1);
291 h = JdkMath.max(h, 1.0e-12 * JdkMath.abs(t0)); // avoids cancellation when computing t1 - t0
292 if (h < getMinStep()) {
293 h = getMinStep();
294 }
295 if (h > getMaxStep()) {
296 h = getMaxStep();
297 }
298 if (! forward) {
299 h = -h;
300 }
301
302 return h;
303 }
304
305 /** Filter the integration step.
306 * @param h signed step
307 * @param forward forward integration indicator
308 * @param acceptSmall if true, steps smaller than the minimal value
309 * are silently increased up to this value, if false such small
310 * steps generate an exception
311 * @return a bounded integration step (h if no bound is reach, or a bounded value)
312 * @exception NumberIsTooSmallException if the step is too small and acceptSmall is false
313 */
314 protected double filterStep(final double h, final boolean forward, final boolean acceptSmall)
315 throws NumberIsTooSmallException {
316
317 double filteredH = h;
318 if (JdkMath.abs(h) < minStep) {
319 if (acceptSmall) {
320 filteredH = forward ? minStep : -minStep;
321 } else {
322 throw new NumberIsTooSmallException(LocalizedFormats.MINIMAL_STEPSIZE_REACHED_DURING_INTEGRATION,
323 JdkMath.abs(h), minStep, true);
324 }
325 }
326
327 if (filteredH > maxStep) {
328 filteredH = maxStep;
329 } else if (filteredH < -maxStep) {
330 filteredH = -maxStep;
331 }
332
333 return filteredH;
334 }
335
336 /** {@inheritDoc} */
337 @Override
338 public abstract void integrate (ExpandableStatefulODE equations, double t)
339 throws NumberIsTooSmallException, DimensionMismatchException,
340 MaxCountExceededException, NoBracketingException;
341
342 /** {@inheritDoc} */
343 @Override
344 public double getCurrentStepStart() {
345 return stepStart;
346 }
347
348 /** Reset internal state to dummy values. */
349 protected void resetInternalState() {
350 stepStart = Double.NaN;
351 stepSize = JdkMath.sqrt(minStep * maxStep);
352 }
353
354 /** Get the minimal step.
355 * @return minimal step
356 */
357 public double getMinStep() {
358 return minStep;
359 }
360
361 /** Get the maximal step.
362 * @return maximal step
363 */
364 public double getMaxStep() {
365 return maxStep;
366 }
367 }