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