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.MathRuntimeException;
21 import org.apache.commons.math.linear.Array2DRowRealMatrix;
22 import org.apache.commons.math.linear.RealMatrix;
23 import org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator;
24 import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator;
25 import org.apache.commons.math.ode.sampling.StepHandler;
26 import org.apache.commons.math.ode.sampling.StepInterpolator;
27
28 /**
29 * This class is the base class for multistep integrators for Ordinary
30 * Differential Equations.
31 * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
32 * <pre>
33 * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
34 * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
35 * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
36 * ...
37 * s<sub>k</sub>(n) = h<sup>k</sup>/k! y(k)<sub>n</sub> for k<sup>th</sup> derivative
38 * </pre></p>
39 * <p>Rather than storing several previous steps separately, this implementation uses
40 * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
41 * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
42 * <pre>
43 * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
44 * </pre>
45 * (we omit the k index in the notation for clarity)</p>
46 * <p>
47 * Multistep integrators with Nordsieck representation are highly sensitive to
48 * large step changes because when the step is multiplied by a factor a, the
49 * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
50 * and the last components are the least accurate ones. The default max growth
51 * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
52 * </p>
53 *
54 * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator
55 * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator
56 * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
57 * @since 2.0
58 */
59 public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
60
61 /** Starter integrator. */
62 private FirstOrderIntegrator starter;
63
64 /** Number of steps of the multistep method (excluding the one being computed). */
65 private final int nSteps;
66
67 /** First scaled derivative (h y'). */
68 protected double[] scaled;
69
70 /** Nordsieck matrix of the higher scaled derivatives.
71 * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y(k))</p>
72 */
73 protected Array2DRowRealMatrix nordsieck;
74
75 /** Stepsize control exponent. */
76 private double exp;
77
78 /** Safety factor for stepsize control. */
79 private double safety;
80
81 /** Minimal reduction factor for stepsize control. */
82 private double minReduction;
83
84 /** Maximal growth factor for stepsize control. */
85 private double maxGrowth;
86
87 /**
88 * Build a multistep integrator with the given stepsize bounds.
89 * <p>The default starter integrator is set to the {@link
90 * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
91 * some defaults settings.</p>
92 * <p>
93 * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
94 * </p>
95 * @param name name of the method
96 * @param nSteps number of steps of the multistep method
97 * (excluding the one being computed)
98 * @param order order of the method
99 * @param minStep minimal step (must be positive even for backward
100 * integration), the last step can be smaller than this
101 * @param maxStep maximal step (must be positive even for backward
102 * integration)
103 * @param scalAbsoluteTolerance allowed absolute error
104 * @param scalRelativeTolerance allowed relative error
105 */
106 protected MultistepIntegrator(final String name, final int nSteps,
107 final int order,
108 final double minStep, final double maxStep,
109 final double scalAbsoluteTolerance,
110 final double scalRelativeTolerance) {
111
112 super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
113
114 if (nSteps <= 0) {
115 throw MathRuntimeException.createIllegalArgumentException(
116 "{0} method needs at least one previous point",
117 name);
118 }
119
120 starter = new DormandPrince853Integrator(minStep, maxStep,
121 scalAbsoluteTolerance,
122 scalRelativeTolerance);
123 this.nSteps = nSteps;
124
125 exp = -1.0 / order;
126
127 // set the default values of the algorithm control parameters
128 setSafety(0.9);
129 setMinReduction(0.2);
130 setMaxGrowth(Math.pow(2.0, -exp));
131
132 }
133
134 /**
135 * Build a multistep integrator with the given stepsize bounds.
136 * <p>The default starter integrator is set to the {@link
137 * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
138 * some defaults settings.</p>
139 * <p>
140 * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
141 * </p>
142 * @param name name of the method
143 * @param nSteps number of steps of the multistep method
144 * (excluding the one being computed)
145 * @param order order of the method
146 * @param minStep minimal step (must be positive even for backward
147 * integration), the last step can be smaller than this
148 * @param maxStep maximal step (must be positive even for backward
149 * integration)
150 * @param vecAbsoluteTolerance allowed absolute error
151 * @param vecRelativeTolerance allowed relative error
152 */
153 protected MultistepIntegrator(final String name, final int nSteps,
154 final int order,
155 final double minStep, final double maxStep,
156 final double[] vecAbsoluteTolerance,
157 final double[] vecRelativeTolerance) {
158 super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
159 starter = new DormandPrince853Integrator(minStep, maxStep,
160 vecAbsoluteTolerance,
161 vecRelativeTolerance);
162 this.nSteps = nSteps;
163
164 exp = -1.0 / order;
165
166 // set the default values of the algorithm control parameters
167 setSafety(0.9);
168 setMinReduction(0.2);
169 setMaxGrowth(Math.pow(2.0, -exp));
170
171 }
172
173 /**
174 * Get the starter integrator.
175 * @return starter integrator
176 */
177 public ODEIntegrator getStarterIntegrator() {
178 return starter;
179 }
180
181 /**
182 * Set the starter integrator.
183 * <p>The various step and event handlers for this starter integrator
184 * will be managed automatically by the multi-step integrator. Any
185 * user configuration for these elements will be cleared before use.</p>
186 * @param starter starter integrator
187 */
188 public void setStarterIntegrator(FirstOrderIntegrator starter) {
189 this.starter = starter;
190 }
191
192 /** Start the integration.
193 * <p>This method computes one step using the underlying starter integrator,
194 * and initializes the Nordsieck vector at step start. The starter integrator
195 * purpose is only to establish initial conditions, it does not really change
196 * time by itself. The top level multistep integrator remains in charge of
197 * handling time propagation and events handling as it will starts its own
198 * computation right from the beginning. In a sense, the starter integrator
199 * can be seen as a dummy one and so it will never trigger any user event nor
200 * call any user step handler.</p>
201 * @param t0 initial time
202 * @param y0 initial value of the state vector at t0
203 * @param t target time for the integration
204 * (can be set to a value smaller than <code>t0</code> for backward integration)
205 * @throws IntegratorException if the integrator cannot perform integration
206 * @throws DerivativeException this exception is propagated to the caller if
207 * the underlying user function triggers one
208 */
209 protected void start(final double t0, final double[] y0, final double t)
210 throws DerivativeException, IntegratorException {
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(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 (DerivativeException de) {
226 if (!(de instanceof InitializationCompletedMarkerException)) {
227 // this is not the expected nominal interruption of the start integrator
228 throw de;
229 }
230 }
231
232 // remove the specific step handler
233 starter.clearStepHandlers();
234
235 }
236
237 /** Initialize the high order scaled derivatives at step start.
238 * @param first first scaled derivative at step start
239 * @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1)
240 * will be modified
241 * @return high order scaled derivatives at step start
242 */
243 protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first,
244 final double[][] multistep);
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 Math.min(maxGrowth, Math.max(minReduction, safety * Math.pow(error, exp)));
294 }
295
296 /** Transformer used to convert the first step to Nordsieck representation. */
297 public static interface NordsieckTransformer {
298 /** Initialize the high order scaled derivatives at step start.
299 * @param first first scaled derivative at step start
300 * @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1)
301 * will be modified
302 * @return high order derivatives at step start
303 */
304 RealMatrix initializeHighOrderDerivatives(double[] first, double[][] multistep);
305 }
306
307 /** Specialized step handler storing the first step. */
308 private class NordsieckInitializer implements StepHandler {
309
310 /** Problem dimension. */
311 private final int n;
312
313 /** Simple constructor.
314 * @param n problem dimension
315 */
316 public NordsieckInitializer(final int n) {
317 this.n = n;
318 }
319
320 /** {@inheritDoc} */
321 public void handleStep(StepInterpolator interpolator, boolean isLast)
322 throws DerivativeException {
323
324 final double prev = interpolator.getPreviousTime();
325 final double curr = interpolator.getCurrentTime();
326 stepStart = prev;
327 stepSize = (curr - prev) / (nSteps + 1);
328
329 // compute the first scaled derivative
330 interpolator.setInterpolatedTime(prev);
331 scaled = interpolator.getInterpolatedDerivatives().clone();
332 for (int j = 0; j < n; ++j) {
333 scaled[j] *= stepSize;
334 }
335
336 // compute the high order scaled derivatives
337 final double[][] multistep = new double[nSteps][];
338 for (int i = 1; i <= nSteps; ++i) {
339 interpolator.setInterpolatedTime(prev + stepSize * i);
340 final double[] msI = interpolator.getInterpolatedDerivatives().clone();
341 for (int j = 0; j < n; ++j) {
342 msI[j] *= stepSize;
343 }
344 multistep[i - 1] = msI;
345 }
346 nordsieck = initializeHighOrderDerivatives(scaled, multistep);
347
348 // stop the integrator after the first step has been handled
349 throw new InitializationCompletedMarkerException();
350
351 }
352
353 /** {@inheritDoc} */
354 public boolean requiresDenseOutput() {
355 return true;
356 }
357
358 /** {@inheritDoc} */
359 public void reset() {
360 // nothing to do
361 }
362
363 }
364
365 /** Marker exception used ONLY to stop the starter integrator after first step. */
366 private static class InitializationCompletedMarkerException
367 extends DerivativeException {
368
369 /** Serializable version identifier. */
370 private static final long serialVersionUID = -4105805787353488365L;
371
372 /** Simple constructor. */
373 public InitializationCompletedMarkerException() {
374 super((Throwable) null);
375 }
376
377 }
378
379 /** Wrapper for differential equations, ensuring start evaluations are counted. */
380 private class CountingDifferentialEquations implements FirstOrderDifferentialEquations {
381
382 /** Dimension of the problem. */
383 private final int dimension;
384
385 /** Simple constructor.
386 * @param dimension dimension of the problem
387 */
388 public CountingDifferentialEquations(final int dimension) {
389 this.dimension = dimension;
390 }
391
392 /** {@inheritDoc} */
393 public void computeDerivatives(double t, double[] y, double[] dot)
394 throws DerivativeException {
395 MultistepIntegrator.this.computeDerivatives(t, y, dot);
396 }
397
398 /** {@inheritDoc} */
399 public int getDimension() {
400 return dimension;
401 }
402 }
403
404 }