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;
019
020 import org.apache.commons.math.MathRuntimeException;
021 import org.apache.commons.math.ode.DerivativeException;
022 import org.apache.commons.math.exception.util.LocalizedFormats;
023 import org.apache.commons.math.linear.Array2DRowRealMatrix;
024 import org.apache.commons.math.linear.RealMatrix;
025 import org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator;
026 import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator;
027 import org.apache.commons.math.ode.sampling.StepHandler;
028 import org.apache.commons.math.ode.sampling.StepInterpolator;
029 import org.apache.commons.math.util.FastMath;
030
031 /**
032 * This class is the base class for multistep integrators for Ordinary
033 * Differential Equations.
034 * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
035 * <pre>
036 * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
037 * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
038 * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
039 * ...
040 * s<sub>k</sub>(n) = h<sup>k</sup>/k! y(k)<sub>n</sub> for k<sup>th</sup> derivative
041 * </pre></p>
042 * <p>Rather than storing several previous steps separately, this implementation uses
043 * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
044 * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
045 * <pre>
046 * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
047 * </pre>
048 * (we omit the k index in the notation for clarity)</p>
049 * <p>
050 * Multistep integrators with Nordsieck representation are highly sensitive to
051 * large step changes because when the step is multiplied by a factor a, the
052 * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
053 * and the last components are the least accurate ones. The default max growth
054 * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
055 * </p>
056 *
057 * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator
058 * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator
059 * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 f??vr. 2011) $
060 * @since 2.0
061 */
062 public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
063
064 /** First scaled derivative (h y'). */
065 protected double[] scaled;
066
067 /** Nordsieck matrix of the higher scaled derivatives.
068 * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y(k))</p>
069 */
070 protected Array2DRowRealMatrix nordsieck;
071
072 /** Starter integrator. */
073 private FirstOrderIntegrator starter;
074
075 /** Number of steps of the multistep method (excluding the one being computed). */
076 private final int nSteps;
077
078 /** Stepsize control exponent. */
079 private double exp;
080
081 /** Safety factor for stepsize control. */
082 private double safety;
083
084 /** Minimal reduction factor for stepsize control. */
085 private double minReduction;
086
087 /** Maximal growth factor for stepsize control. */
088 private double maxGrowth;
089
090 /**
091 * Build a multistep integrator with the given stepsize bounds.
092 * <p>The default starter integrator is set to the {@link
093 * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
094 * some defaults settings.</p>
095 * <p>
096 * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
097 * </p>
098 * @param name name of the method
099 * @param nSteps number of steps of the multistep method
100 * (excluding the one being computed)
101 * @param order order of the method
102 * @param minStep minimal step (must be positive even for backward
103 * integration), the last step can be smaller than this
104 * @param maxStep maximal step (must be positive even for backward
105 * integration)
106 * @param scalAbsoluteTolerance allowed absolute error
107 * @param scalRelativeTolerance allowed relative error
108 */
109 protected MultistepIntegrator(final String name, final int nSteps,
110 final int order,
111 final double minStep, final double maxStep,
112 final double scalAbsoluteTolerance,
113 final double scalRelativeTolerance) {
114
115 super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
116
117 if (nSteps <= 0) {
118 throw MathRuntimeException.createIllegalArgumentException(
119 LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_ONE_PREVIOUS_POINT,
120 name);
121 }
122
123 starter = new DormandPrince853Integrator(minStep, maxStep,
124 scalAbsoluteTolerance,
125 scalRelativeTolerance);
126 this.nSteps = nSteps;
127
128 exp = -1.0 / order;
129
130 // set the default values of the algorithm control parameters
131 setSafety(0.9);
132 setMinReduction(0.2);
133 setMaxGrowth(FastMath.pow(2.0, -exp));
134
135 }
136
137 /**
138 * Build a multistep integrator with the given stepsize bounds.
139 * <p>The default starter integrator is set to the {@link
140 * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
141 * some defaults settings.</p>
142 * <p>
143 * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
144 * </p>
145 * @param name name of the method
146 * @param nSteps number of steps of the multistep method
147 * (excluding the one being computed)
148 * @param order order of the method
149 * @param minStep minimal step (must be positive even for backward
150 * integration), the last step can be smaller than this
151 * @param maxStep maximal step (must be positive even for backward
152 * integration)
153 * @param vecAbsoluteTolerance allowed absolute error
154 * @param vecRelativeTolerance allowed relative error
155 */
156 protected MultistepIntegrator(final String name, final int nSteps,
157 final int order,
158 final double minStep, final double maxStep,
159 final double[] vecAbsoluteTolerance,
160 final double[] vecRelativeTolerance) {
161 super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
162 starter = new DormandPrince853Integrator(minStep, maxStep,
163 vecAbsoluteTolerance,
164 vecRelativeTolerance);
165 this.nSteps = nSteps;
166
167 exp = -1.0 / order;
168
169 // set the default values of the algorithm control parameters
170 setSafety(0.9);
171 setMinReduction(0.2);
172 setMaxGrowth(FastMath.pow(2.0, -exp));
173
174 }
175
176 /**
177 * Get the starter integrator.
178 * @return starter integrator
179 */
180 public ODEIntegrator getStarterIntegrator() {
181 return starter;
182 }
183
184 /**
185 * Set the starter integrator.
186 * <p>The various step and event handlers for this starter integrator
187 * will be managed automatically by the multi-step integrator. Any
188 * user configuration for these elements will be cleared before use.</p>
189 * @param starterIntegrator starter integrator
190 */
191 public void setStarterIntegrator(FirstOrderIntegrator starterIntegrator) {
192 this.starter = starterIntegrator;
193 }
194
195 /** Start the integration.
196 * <p>This method computes one step using the underlying starter integrator,
197 * and initializes the Nordsieck vector at step start. The starter integrator
198 * purpose is only to establish initial conditions, it does not really change
199 * time by itself. The top level multistep integrator remains in charge of
200 * handling time propagation and events handling as it will starts its own
201 * computation right from the beginning. In a sense, the starter integrator
202 * can be seen as a dummy one and so it will never trigger any user event nor
203 * call any user step handler.</p>
204 * @param t0 initial time
205 * @param y0 initial value of the state vector at t0
206 * @param t target time for the integration
207 * (can be set to a value smaller than <code>t0</code> for backward integration)
208 * @throws IntegratorException if the integrator cannot perform integration
209 * @throws DerivativeException this exception is propagated to the caller if
210 * the underlying user function triggers one
211 */
212 protected void start(final double t0, final double[] y0, final double t)
213 throws DerivativeException, IntegratorException {
214
215 // make sure NO user event nor user step handler is triggered,
216 // this is the task of the top level integrator, not the task
217 // of the starter integrator
218 starter.clearEventHandlers();
219 starter.clearStepHandlers();
220
221 // set up one specific step handler to extract initial Nordsieck vector
222 starter.addStepHandler(new NordsieckInitializer(y0.length));
223
224 // start integration, expecting a InitializationCompletedMarkerException
225 try {
226 starter.integrate(new CountingDifferentialEquations(y0.length),
227 t0, y0, t, new double[y0.length]);
228 } catch (DerivativeException mue) {
229 if (!(mue instanceof InitializationCompletedMarkerException)) {
230 // this is not the expected nominal interruption of the start integrator
231 throw mue;
232 }
233 }
234
235 // remove the specific step handler
236 starter.clearStepHandlers();
237
238 }
239
240 /** Initialize the high order scaled derivatives at step start.
241 * @param first first scaled derivative at step start
242 * @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1)
243 * will be modified
244 * @return high order scaled derivatives at step start
245 */
246 protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first,
247 final double[][] multistep);
248
249 /** Get the minimal reduction factor for stepsize control.
250 * @return minimal reduction factor
251 */
252 public double getMinReduction() {
253 return minReduction;
254 }
255
256 /** Set the minimal reduction factor for stepsize control.
257 * @param minReduction minimal reduction factor
258 */
259 public void setMinReduction(final double minReduction) {
260 this.minReduction = minReduction;
261 }
262
263 /** Get the maximal growth factor for stepsize control.
264 * @return maximal growth factor
265 */
266 public double getMaxGrowth() {
267 return maxGrowth;
268 }
269
270 /** Set the maximal growth factor for stepsize control.
271 * @param maxGrowth maximal growth factor
272 */
273 public void setMaxGrowth(final double maxGrowth) {
274 this.maxGrowth = maxGrowth;
275 }
276
277 /** Get the safety factor for stepsize control.
278 * @return safety factor
279 */
280 public double getSafety() {
281 return safety;
282 }
283
284 /** Set the safety factor for stepsize control.
285 * @param safety safety factor
286 */
287 public void setSafety(final double safety) {
288 this.safety = safety;
289 }
290
291 /** Compute step grow/shrink factor according to normalized error.
292 * @param error normalized error of the current step
293 * @return grow/shrink factor for next step
294 */
295 protected double computeStepGrowShrinkFactor(final double error) {
296 return FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
297 }
298
299 /** Transformer used to convert the first step to Nordsieck representation. */
300 public static interface NordsieckTransformer {
301 /** Initialize the high order scaled derivatives at step start.
302 * @param first first scaled derivative at step start
303 * @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1)
304 * will be modified
305 * @return high order derivatives at step start
306 */
307 RealMatrix initializeHighOrderDerivatives(double[] first, double[][] multistep);
308 }
309
310 /** Specialized step handler storing the first step. */
311 private class NordsieckInitializer implements StepHandler {
312
313 /** Problem dimension. */
314 private final int n;
315
316 /** Simple constructor.
317 * @param n problem dimension
318 */
319 public NordsieckInitializer(final int n) {
320 this.n = n;
321 }
322
323 /** {@inheritDoc} */
324 public void handleStep(StepInterpolator interpolator, boolean isLast)
325 throws DerivativeException {
326
327 final double prev = interpolator.getPreviousTime();
328 final double curr = interpolator.getCurrentTime();
329 stepStart = prev;
330 stepSize = (curr - prev) / (nSteps + 1);
331
332 // compute the first scaled derivative
333 interpolator.setInterpolatedTime(prev);
334 scaled = interpolator.getInterpolatedDerivatives().clone();
335 for (int j = 0; j < n; ++j) {
336 scaled[j] *= stepSize;
337 }
338
339 // compute the high order scaled derivatives
340 final double[][] multistep = new double[nSteps][];
341 for (int i = 1; i <= nSteps; ++i) {
342 interpolator.setInterpolatedTime(prev + stepSize * i);
343 final double[] msI = interpolator.getInterpolatedDerivatives().clone();
344 for (int j = 0; j < n; ++j) {
345 msI[j] *= stepSize;
346 }
347 multistep[i - 1] = msI;
348 }
349 nordsieck = initializeHighOrderDerivatives(scaled, multistep);
350
351 // stop the integrator after the first step has been handled
352 throw new InitializationCompletedMarkerException();
353
354 }
355
356 /** {@inheritDoc} */
357 public boolean requiresDenseOutput() {
358 return true;
359 }
360
361 /** {@inheritDoc} */
362 public void reset() {
363 // nothing to do
364 }
365
366 }
367
368 /** Marker exception used ONLY to stop the starter integrator after first step. */
369 private static class InitializationCompletedMarkerException
370 extends DerivativeException {
371
372 /** Serializable version identifier. */
373 private static final long serialVersionUID = -4105805787353488365L;
374
375 /** Simple constructor. */
376 public InitializationCompletedMarkerException() {
377 super((Throwable) null);
378 }
379
380 }
381
382 /** Wrapper for differential equations, ensuring start evaluations are counted. */
383 private class CountingDifferentialEquations implements ExtendedFirstOrderDifferentialEquations {
384
385 /** Dimension of the problem. */
386 private final int dimension;
387
388 /** Simple constructor.
389 * @param dimension dimension of the problem
390 */
391 public CountingDifferentialEquations(final int dimension) {
392 this.dimension = dimension;
393 }
394
395 /** {@inheritDoc} */
396 public void computeDerivatives(double t, double[] y, double[] dot)
397 throws DerivativeException {
398 MultistepIntegrator.this.computeDerivatives(t, y, dot);
399 }
400
401 /** {@inheritDoc} */
402 public int getDimension() {
403 return dimension;
404 }
405
406 /** {@inheritDoc} */
407 public int getMainSetDimension() {
408 return mainSetDimension;
409 }
410 }
411
412 }