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