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.math.ode;
19
20 import org.apache.commons.math.exception.MathIllegalArgumentException;
21 import org.apache.commons.math.exception.MathIllegalStateException;
22 import org.apache.commons.math.exception.util.LocalizedFormats;
23 import org.apache.commons.math.linear.Array2DRowRealMatrix;
24 import org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator;
25 import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator;
26 import org.apache.commons.math.ode.sampling.StepHandler;
27 import org.apache.commons.math.ode.sampling.StepInterpolator;
28 import org.apache.commons.math.util.FastMath;
29
30 /**
31 * This class is the base class for multistep integrators for Ordinary
32 * Differential Equations.
33 * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
34 * <pre>
35 * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
36 * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
37 * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
38 * ...
39 * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
40 * </pre></p>
41 * <p>Rather than storing several previous steps separately, this implementation uses
42 * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
43 * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
44 * <pre>
45 * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
46 * </pre>
47 * (we omit the k index in the notation for clarity)</p>
48 * <p>
49 * Multistep integrators with Nordsieck representation are highly sensitive to
50 * large step changes because when the step is multiplied by factor a, the
51 * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
52 * and the last components are the least accurate ones. The default max growth
53 * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
54 * </p>
55 *
56 * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator
57 * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator
58 * @version $Id: MultistepIntegrator.java 1178235 2011-10-02 19:43:17Z luc $
59 * @since 2.0
60 */
61 public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
62
63 /** First scaled derivative (h y'). */
64 protected double[] scaled;
65
66 /** Nordsieck matrix of the higher scaled derivatives.
67 * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y<sup>(k)</sup>)</p>
68 */
69 protected Array2DRowRealMatrix nordsieck;
70
71 /** Starter integrator. */
72 private FirstOrderIntegrator starter;
73
74 /** Number of steps of the multistep method (excluding the one being computed). */
75 private final int nSteps;
76
77 /** Stepsize control exponent. */
78 private double exp;
79
80 /** Safety factor for stepsize control. */
81 private double safety;
82
83 /** Minimal reduction factor for stepsize control. */
84 private double minReduction;
85
86 /** Maximal growth factor for stepsize control. */
87 private double maxGrowth;
88
89 /**
90 * Build a multistep integrator with the given stepsize bounds.
91 * <p>The default starter integrator is set to the {@link
92 * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
93 * some defaults settings.</p>
94 * <p>
95 * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
96 * </p>
97 * @param name name of the method
98 * @param nSteps number of steps of the multistep method
99 * (excluding the one being computed)
100 * @param order order of the method
101 * @param minStep minimal step (must be positive even for backward
102 * integration), the last step can be smaller than this
103 * @param maxStep maximal step (must be positive even for backward
104 * integration)
105 * @param scalAbsoluteTolerance allowed absolute error
106 * @param scalRelativeTolerance allowed relative error
107 */
108 protected MultistepIntegrator(final String name, final int nSteps,
109 final int order,
110 final double minStep, final double maxStep,
111 final double scalAbsoluteTolerance,
112 final double scalRelativeTolerance) {
113
114 super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
115
116 if (nSteps <= 1) {
117 throw new MathIllegalArgumentException(
118 LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
119 name);
120 }
121
122 starter = new DormandPrince853Integrator(minStep, maxStep,
123 scalAbsoluteTolerance,
124 scalRelativeTolerance);
125 this.nSteps = nSteps;
126
127 exp = -1.0 / order;
128
129 // set the default values of the algorithm control parameters
130 setSafety(0.9);
131 setMinReduction(0.2);
132 setMaxGrowth(FastMath.pow(2.0, -exp));
133
134 }
135
136 /**
137 * Build a multistep integrator with the given stepsize bounds.
138 * <p>The default starter integrator is set to the {@link
139 * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
140 * some defaults settings.</p>
141 * <p>
142 * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
143 * </p>
144 * @param name name of the method
145 * @param nSteps number of steps of the multistep method
146 * (excluding the one being computed)
147 * @param order order of the method
148 * @param minStep minimal step (must be positive even for backward
149 * integration), the last step can be smaller than this
150 * @param maxStep maximal step (must be positive even for backward
151 * integration)
152 * @param vecAbsoluteTolerance allowed absolute error
153 * @param vecRelativeTolerance allowed relative error
154 */
155 protected MultistepIntegrator(final String name, final int nSteps,
156 final int order,
157 final double minStep, final double maxStep,
158 final double[] vecAbsoluteTolerance,
159 final double[] vecRelativeTolerance) {
160 super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
161 starter = new DormandPrince853Integrator(minStep, maxStep,
162 vecAbsoluteTolerance,
163 vecRelativeTolerance);
164 this.nSteps = nSteps;
165
166 exp = -1.0 / order;
167
168 // set the default values of the algorithm control parameters
169 setSafety(0.9);
170 setMinReduction(0.2);
171 setMaxGrowth(FastMath.pow(2.0, -exp));
172
173 }
174
175 /**
176 * Get the starter integrator.
177 * @return starter integrator
178 */
179 public ODEIntegrator getStarterIntegrator() {
180 return starter;
181 }
182
183 /**
184 * Set the starter integrator.
185 * <p>The various step and event handlers for this starter integrator
186 * will be managed automatically by the multi-step integrator. Any
187 * user configuration for these elements will be cleared before use.</p>
188 * @param starterIntegrator starter integrator
189 */
190 public void setStarterIntegrator(FirstOrderIntegrator starterIntegrator) {
191 this.starter = starterIntegrator;
192 }
193
194 /** Start the integration.
195 * <p>This method computes one step using the underlying starter integrator,
196 * and initializes the Nordsieck vector at step start. The starter integrator
197 * purpose is only to establish initial conditions, it does not really change
198 * time by itself. The top level multistep integrator remains in charge of
199 * handling time propagation and events handling as it will starts its own
200 * computation right from the beginning. In a sense, the starter integrator
201 * can be seen as a dummy one and so it will never trigger any user event nor
202 * call any user step handler.</p>
203 * @param t0 initial time
204 * @param y0 initial value of the state vector at t0
205 * @param t target time for the integration
206 * (can be set to a value smaller than <code>t0</code> for backward integration)
207 * @throws MathIllegalStateException if the integrator cannot perform integration
208 */
209 protected void start(final double t0, final double[] y0, final double t)
210 throws MathIllegalStateException {
211
212 // make sure NO user event nor user step handler is triggered,
213 // this is the task of the top level integrator, not the task
214 // of the starter integrator
215 starter.clearEventHandlers();
216 starter.clearStepHandlers();
217
218 // set up one specific step handler to extract initial Nordsieck vector
219 starter.addStepHandler(new NordsieckInitializer(nSteps, y0.length));
220
221 // start integration, expecting a InitializationCompletedMarkerException
222 try {
223 starter.integrate(new CountingDifferentialEquations(y0.length),
224 t0, y0, t, new double[y0.length]);
225 } catch (InitializationCompletedMarkerException icme) {
226 // this is the expected nominal interruption of the start integrator
227 }
228
229 // remove the specific step handler
230 starter.clearStepHandlers();
231
232 }
233
234 /** Initialize the high order scaled derivatives at step start.
235 * @param h step size to use for scaling
236 * @param t first steps times
237 * @param y first steps states
238 * @param yDot first steps derivatives
239 * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
240 * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
241 */
242 protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
243 final double[][] y,
244 final double[][] yDot);
245
246 /** Get the minimal reduction factor for stepsize control.
247 * @return minimal reduction factor
248 */
249 public double getMinReduction() {
250 return minReduction;
251 }
252
253 /** Set the minimal reduction factor for stepsize control.
254 * @param minReduction minimal reduction factor
255 */
256 public void setMinReduction(final double minReduction) {
257 this.minReduction = minReduction;
258 }
259
260 /** Get the maximal growth factor for stepsize control.
261 * @return maximal growth factor
262 */
263 public double getMaxGrowth() {
264 return maxGrowth;
265 }
266
267 /** Set the maximal growth factor for stepsize control.
268 * @param maxGrowth maximal growth factor
269 */
270 public void setMaxGrowth(final double maxGrowth) {
271 this.maxGrowth = maxGrowth;
272 }
273
274 /** Get the safety factor for stepsize control.
275 * @return safety factor
276 */
277 public double getSafety() {
278 return safety;
279 }
280
281 /** Set the safety factor for stepsize control.
282 * @param safety safety factor
283 */
284 public void setSafety(final double safety) {
285 this.safety = safety;
286 }
287
288 /** Compute step grow/shrink factor according to normalized error.
289 * @param error normalized error of the current step
290 * @return grow/shrink factor for next step
291 */
292 protected double computeStepGrowShrinkFactor(final double error) {
293 return FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
294 }
295
296 /** Transformer used to convert the first step to Nordsieck representation. */
297 public interface NordsieckTransformer {
298 /** Initialize the high order scaled derivatives at step start.
299 * @param h step size to use for scaling
300 * @param t first steps times
301 * @param y first steps states
302 * @param yDot first steps derivatives
303 * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
304 * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
305 */
306 Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
307 final double[][] y,
308 final double[][] yDot);
309 }
310
311 /** Specialized step handler storing the first step. */
312 private class NordsieckInitializer implements StepHandler {
313
314 /** Steps counter. */
315 private int count;
316
317 /** First steps times. */
318 private final double[] t;
319
320 /** First steps states. */
321 private final double[][] y;
322
323 /** First steps derivatives. */
324 private final double[][] yDot;
325
326 /** Simple constructor.
327 * @param nSteps number of steps of the multistep method (excluding the one being computed)
328 * @param n problem dimension
329 */
330 public NordsieckInitializer(final int nSteps, final int n) {
331 this.count = 0;
332 this.t = new double[nSteps];
333 this.y = new double[nSteps][n];
334 this.yDot = new double[nSteps][n];
335 }
336
337 /** {@inheritDoc} */
338 public void handleStep(StepInterpolator interpolator, boolean isLast) {
339
340 final double prev = interpolator.getPreviousTime();
341 final double curr = interpolator.getCurrentTime();
342
343 if (count == 0) {
344 // first step, we need to store also the beginning of the step
345 interpolator.setInterpolatedTime(prev);
346 t[0] = prev;
347 System.arraycopy(interpolator.getInterpolatedState(), 0,
348 y[0], 0, y[0].length);
349 System.arraycopy(interpolator.getInterpolatedDerivatives(), 0,
350 yDot[0], 0, yDot[0].length);
351 }
352
353 // store the end of the step
354 ++count;
355 interpolator.setInterpolatedTime(curr);
356 t[count] = curr;
357 System.arraycopy(interpolator.getInterpolatedState(), 0,
358 y[count], 0, y[count].length);
359 System.arraycopy(interpolator.getInterpolatedDerivatives(), 0,
360 yDot[count], 0, yDot[count].length);
361
362 if (count == t.length - 1) {
363
364 // this was the last step we needed, we can compute the derivatives
365 stepStart = t[0];
366 stepSize = (t[t.length - 1] - t[0]) / (t.length - 1);
367
368 // first scaled derivative
369 scaled = yDot[0].clone();
370 for (int j = 0; j < scaled.length; ++j) {
371 scaled[j] *= stepSize;
372 }
373
374 // higher order derivatives
375 nordsieck = initializeHighOrderDerivatives(stepSize, t, y, yDot);
376
377 // stop the integrator now that all needed steps have been handled
378 throw new InitializationCompletedMarkerException();
379
380 }
381
382 }
383
384 /** {@inheritDoc} */
385 public void reset() {
386 // nothing to do
387 }
388
389 }
390
391 /** Marker exception used ONLY to stop the starter integrator after first step. */
392 private static class InitializationCompletedMarkerException
393 extends RuntimeException {
394
395 /** Serializable version identifier. */
396 private static final long serialVersionUID = -1914085471038046418L;
397
398 /** Simple constructor. */
399 public InitializationCompletedMarkerException() {
400 super((Throwable) null);
401 }
402
403 }
404
405 /** Wrapper for differential equations, ensuring start evaluations are counted. */
406 private class CountingDifferentialEquations implements FirstOrderDifferentialEquations {
407
408 /** Dimension of the problem. */
409 private final int dimension;
410
411 /** Simple constructor.
412 * @param dimension dimension of the problem
413 */
414 public CountingDifferentialEquations(final int dimension) {
415 this.dimension = dimension;
416 }
417
418 /** {@inheritDoc} */
419 public void computeDerivatives(double t, double[] y, double[] dot) {
420 MultistepIntegrator.this.computeDerivatives(t, y, dot);
421 }
422
423 /** {@inheritDoc} */
424 public int getDimension() {
425 return dimension;
426 }
427
428 }
429
430 }