| /* |
| * Licensed to the Apache Software Foundation (ASF) under one or more |
| * contributor license agreements. See the NOTICE file distributed with |
| * this work for additional information regarding copyright ownership. |
| * The ASF licenses this file to You under the Apache License, Version 2.0 |
| * (the "License"); you may not use this file except in compliance with |
| * the License. You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| |
| package org.apache.commons.math4.legacy.ode; |
| |
| import org.apache.commons.math4.legacy.core.Field; |
| import org.apache.commons.math4.legacy.core.RealFieldElement; |
| import org.apache.commons.math4.legacy.exception.DimensionMismatchException; |
| import org.apache.commons.math4.legacy.exception.MathIllegalStateException; |
| import org.apache.commons.math4.legacy.exception.MaxCountExceededException; |
| import org.apache.commons.math4.legacy.exception.NoBracketingException; |
| import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException; |
| import org.apache.commons.math4.legacy.exception.util.LocalizedFormats; |
| import org.apache.commons.math4.legacy.linear.Array2DRowFieldMatrix; |
| import org.apache.commons.math4.legacy.ode.nonstiff.AdaptiveStepsizeFieldIntegrator; |
| import org.apache.commons.math4.legacy.ode.nonstiff.DormandPrince853FieldIntegrator; |
| import org.apache.commons.math4.legacy.ode.sampling.FieldStepHandler; |
| import org.apache.commons.math4.legacy.ode.sampling.FieldStepInterpolator; |
| import org.apache.commons.math4.legacy.core.jdkmath.AccurateMath; |
| import org.apache.commons.math4.legacy.core.MathArrays; |
| |
| /** |
| * This class is the base class for multistep integrators for Ordinary |
| * Differential Equations. |
| * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as: |
| * <div style="white-space: pre"><code> |
| * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative |
| * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative |
| * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative |
| * ... |
| * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative |
| * </code></div> |
| * <p>Rather than storing several previous steps separately, this implementation uses |
| * the Nordsieck vector with higher degrees scaled derivatives all taken at the same |
| * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as: |
| * <div style="white-space: pre"><code> |
| * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup> |
| * </code></div> |
| * (we omit the k index in the notation for clarity) |
| * <p> |
| * Multistep integrators with Nordsieck representation are highly sensitive to |
| * large step changes because when the step is multiplied by factor a, the |
| * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup> |
| * and the last components are the least accurate ones. The default max growth |
| * factor is therefore set to a quite low value: 2<sup>1/order</sup>. |
| * </p> |
| * |
| * @see org.apache.commons.math4.legacy.ode.nonstiff.AdamsBashforthFieldIntegrator |
| * @see org.apache.commons.math4.legacy.ode.nonstiff.AdamsMoultonFieldIntegrator |
| * @param <T> the type of the field elements |
| * @since 3.6 |
| */ |
| public abstract class MultistepFieldIntegrator<T extends RealFieldElement<T>> |
| extends AdaptiveStepsizeFieldIntegrator<T> { |
| |
| /** First scaled derivative (h y'). */ |
| protected T[] scaled; |
| |
| /** Nordsieck matrix of the higher scaled derivatives. |
| * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y<sup>(k)</sup>)</p> |
| */ |
| protected Array2DRowFieldMatrix<T> nordsieck; |
| |
| /** Starter integrator. */ |
| private FirstOrderFieldIntegrator<T> starter; |
| |
| /** Number of steps of the multistep method (excluding the one being computed). */ |
| private final int nSteps; |
| |
| /** Stepsize control exponent. */ |
| private double exp; |
| |
| /** Safety factor for stepsize control. */ |
| private double safety; |
| |
| /** Minimal reduction factor for stepsize control. */ |
| private double minReduction; |
| |
| /** Maximal growth factor for stepsize control. */ |
| private double maxGrowth; |
| |
| /** |
| * Build a multistep integrator with the given stepsize bounds. |
| * <p>The default starter integrator is set to the {@link |
| * DormandPrince853FieldIntegrator Dormand-Prince 8(5,3)} integrator with |
| * some defaults settings.</p> |
| * <p> |
| * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>. |
| * </p> |
| * @param field field to which the time and state vector elements belong |
| * @param name name of the method |
| * @param nSteps number of steps of the multistep method |
| * (excluding the one being computed) |
| * @param order order of the method |
| * @param minStep minimal step (must be positive even for backward |
| * integration), the last step can be smaller than this |
| * @param maxStep maximal step (must be positive even for backward |
| * integration) |
| * @param scalAbsoluteTolerance allowed absolute error |
| * @param scalRelativeTolerance allowed relative error |
| * @exception NumberIsTooSmallException if number of steps is smaller than 2 |
| */ |
| protected MultistepFieldIntegrator(final Field<T> field, final String name, |
| final int nSteps, final int order, |
| final double minStep, final double maxStep, |
| final double scalAbsoluteTolerance, |
| final double scalRelativeTolerance) |
| throws NumberIsTooSmallException { |
| |
| super(field, name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); |
| |
| if (nSteps < 2) { |
| throw new NumberIsTooSmallException( |
| LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS, |
| nSteps, 2, true); |
| } |
| |
| starter = new DormandPrince853FieldIntegrator<>(field, minStep, maxStep, |
| scalAbsoluteTolerance, |
| scalRelativeTolerance); |
| this.nSteps = nSteps; |
| |
| exp = -1.0 / order; |
| |
| // set the default values of the algorithm control parameters |
| setSafety(0.9); |
| setMinReduction(0.2); |
| setMaxGrowth(AccurateMath.pow(2.0, -exp)); |
| |
| } |
| |
| /** |
| * Build a multistep integrator with the given stepsize bounds. |
| * <p>The default starter integrator is set to the {@link |
| * DormandPrince853FieldIntegrator Dormand-Prince 8(5,3)} integrator with |
| * some defaults settings.</p> |
| * <p> |
| * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>. |
| * </p> |
| * @param field field to which the time and state vector elements belong |
| * @param name name of the method |
| * @param nSteps number of steps of the multistep method |
| * (excluding the one being computed) |
| * @param order order of the method |
| * @param minStep minimal step (must be positive even for backward |
| * integration), the last step can be smaller than this |
| * @param maxStep maximal step (must be positive even for backward |
| * integration) |
| * @param vecAbsoluteTolerance allowed absolute error |
| * @param vecRelativeTolerance allowed relative error |
| */ |
| protected MultistepFieldIntegrator(final Field<T> field, final String name, final int nSteps, |
| final int order, |
| final double minStep, final double maxStep, |
| final double[] vecAbsoluteTolerance, |
| final double[] vecRelativeTolerance) { |
| super(field, name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); |
| starter = new DormandPrince853FieldIntegrator<>(field, minStep, maxStep, |
| vecAbsoluteTolerance, |
| vecRelativeTolerance); |
| this.nSteps = nSteps; |
| |
| exp = -1.0 / order; |
| |
| // set the default values of the algorithm control parameters |
| setSafety(0.9); |
| setMinReduction(0.2); |
| setMaxGrowth(AccurateMath.pow(2.0, -exp)); |
| |
| } |
| |
| /** |
| * Get the starter integrator. |
| * @return starter integrator |
| */ |
| public FirstOrderFieldIntegrator<T> getStarterIntegrator() { |
| return starter; |
| } |
| |
| /** |
| * Set the starter integrator. |
| * <p>The various step and event handlers for this starter integrator |
| * will be managed automatically by the multi-step integrator. Any |
| * user configuration for these elements will be cleared before use.</p> |
| * @param starterIntegrator starter integrator |
| */ |
| public void setStarterIntegrator(FirstOrderFieldIntegrator<T> starterIntegrator) { |
| this.starter = starterIntegrator; |
| } |
| |
| /** Start the integration. |
| * <p>This method computes one step using the underlying starter integrator, |
| * and initializes the Nordsieck vector at step start. The starter integrator |
| * purpose is only to establish initial conditions, it does not really change |
| * time by itself. The top level multistep integrator remains in charge of |
| * handling time propagation and events handling as it will starts its own |
| * computation right from the beginning. In a sense, the starter integrator |
| * can be seen as a dummy one and so it will never trigger any user event nor |
| * call any user step handler.</p> |
| * @param equations complete set of differential equations to integrate |
| * @param initialState initial state (time, primary and secondary state vectors) |
| * @param t target time for the integration |
| * (can be set to a value smaller than <code>t0</code> for backward integration) |
| * @exception DimensionMismatchException if arrays dimension do not match equations settings |
| * @exception NumberIsTooSmallException if integration step is too small |
| * @exception MaxCountExceededException if the number of functions evaluations is exceeded |
| * @exception NoBracketingException if the location of an event cannot be bracketed |
| */ |
| protected void start(final FieldExpandableODE<T> equations, final FieldODEState<T> initialState, final T t) |
| throws DimensionMismatchException, NumberIsTooSmallException, |
| MaxCountExceededException, NoBracketingException { |
| |
| // make sure NO user event nor user step handler is triggered, |
| // this is the task of the top level integrator, not the task |
| // of the starter integrator |
| starter.clearEventHandlers(); |
| starter.clearStepHandlers(); |
| |
| // set up one specific step handler to extract initial Nordsieck vector |
| starter.addStepHandler(new FieldNordsieckInitializer(equations.getMapper(), (nSteps + 3) / 2)); |
| |
| // start integration, expecting a InitializationCompletedMarkerException |
| try { |
| |
| starter.integrate(equations, initialState, t); |
| |
| // we should not reach this step |
| throw new MathIllegalStateException(LocalizedFormats.MULTISTEP_STARTER_STOPPED_EARLY); |
| |
| } catch (InitializationCompletedMarkerException icme) { // NOPMD |
| // this is the expected nominal interruption of the start integrator |
| |
| // count the evaluations used by the starter |
| getEvaluationsCounter().increment(starter.getEvaluations()); |
| |
| } |
| |
| // remove the specific step handler |
| starter.clearStepHandlers(); |
| |
| } |
| |
| /** Initialize the high order scaled derivatives at step start. |
| * @param h step size to use for scaling |
| * @param t first steps times |
| * @param y first steps states |
| * @param yDot first steps derivatives |
| * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>, |
| * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>) |
| */ |
| protected abstract Array2DRowFieldMatrix<T> initializeHighOrderDerivatives(T h, T[] t, |
| T[][] y, |
| T[][] yDot); |
| |
| /** Get the minimal reduction factor for stepsize control. |
| * @return minimal reduction factor |
| */ |
| public double getMinReduction() { |
| return minReduction; |
| } |
| |
| /** Set the minimal reduction factor for stepsize control. |
| * @param minReduction minimal reduction factor |
| */ |
| public void setMinReduction(final double minReduction) { |
| this.minReduction = minReduction; |
| } |
| |
| /** Get the maximal growth factor for stepsize control. |
| * @return maximal growth factor |
| */ |
| public double getMaxGrowth() { |
| return maxGrowth; |
| } |
| |
| /** Set the maximal growth factor for stepsize control. |
| * @param maxGrowth maximal growth factor |
| */ |
| public void setMaxGrowth(final double maxGrowth) { |
| this.maxGrowth = maxGrowth; |
| } |
| |
| /** Get the safety factor for stepsize control. |
| * @return safety factor |
| */ |
| public double getSafety() { |
| return safety; |
| } |
| |
| /** Set the safety factor for stepsize control. |
| * @param safety safety factor |
| */ |
| public void setSafety(final double safety) { |
| this.safety = safety; |
| } |
| |
| /** Get the number of steps of the multistep method (excluding the one being computed). |
| * @return number of steps of the multistep method (excluding the one being computed) |
| */ |
| public int getNSteps() { |
| return nSteps; |
| } |
| |
| /** Rescale the instance. |
| * <p>Since the scaled and Nordsieck arrays are shared with the caller, |
| * this method has the side effect of rescaling this arrays in the caller too.</p> |
| * @param newStepSize new step size to use in the scaled and Nordsieck arrays |
| */ |
| protected void rescale(final T newStepSize) { |
| |
| final T ratio = newStepSize.divide(getStepSize()); |
| for (int i = 0; i < scaled.length; ++i) { |
| scaled[i] = scaled[i].multiply(ratio); |
| } |
| |
| final T[][] nData = nordsieck.getDataRef(); |
| T power = ratio; |
| for (int i = 0; i < nData.length; ++i) { |
| power = power.multiply(ratio); |
| final T[] nDataI = nData[i]; |
| for (int j = 0; j < nDataI.length; ++j) { |
| nDataI[j] = nDataI[j].multiply(power); |
| } |
| } |
| |
| setStepSize(newStepSize); |
| |
| } |
| |
| |
| /** Compute step grow/shrink factor according to normalized error. |
| * @param error normalized error of the current step |
| * @return grow/shrink factor for next step |
| */ |
| protected T computeStepGrowShrinkFactor(final T error) { |
| return RealFieldElement.min(error.getField().getZero().add(maxGrowth), |
| RealFieldElement.max(error.getField().getZero().add(minReduction), |
| error.pow(exp).multiply(safety))); |
| } |
| |
| /** Specialized step handler storing the first step. |
| */ |
| private class FieldNordsieckInitializer implements FieldStepHandler<T> { |
| |
| /** Equation mapper. */ |
| private final FieldEquationsMapper<T> mapper; |
| |
| /** Steps counter. */ |
| private int count; |
| |
| /** Saved start. */ |
| private FieldODEStateAndDerivative<T> savedStart; |
| |
| /** First steps times. */ |
| private final T[] t; |
| |
| /** First steps states. */ |
| private final T[][] y; |
| |
| /** First steps derivatives. */ |
| private final T[][] yDot; |
| |
| /** Simple constructor. |
| * @param mapper equation mapper |
| * @param nbStartPoints number of start points (including the initial point) |
| */ |
| FieldNordsieckInitializer(final FieldEquationsMapper<T> mapper, final int nbStartPoints) { |
| this.mapper = mapper; |
| this.count = 0; |
| this.t = MathArrays.buildArray(getField(), nbStartPoints); |
| this.y = MathArrays.buildArray(getField(), nbStartPoints, -1); |
| this.yDot = MathArrays.buildArray(getField(), nbStartPoints, -1); |
| } |
| |
| /** {@inheritDoc} */ |
| @Override |
| public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast) |
| throws MaxCountExceededException { |
| |
| |
| if (count == 0) { |
| // first step, we need to store also the point at the beginning of the step |
| final FieldODEStateAndDerivative<T> prev = interpolator.getPreviousState(); |
| savedStart = prev; |
| t[count] = prev.getTime(); |
| y[count] = mapper.mapState(prev); |
| yDot[count] = mapper.mapDerivative(prev); |
| } |
| |
| // store the point at the end of the step |
| ++count; |
| final FieldODEStateAndDerivative<T> curr = interpolator.getCurrentState(); |
| t[count] = curr.getTime(); |
| y[count] = mapper.mapState(curr); |
| yDot[count] = mapper.mapDerivative(curr); |
| |
| if (count == t.length - 1) { |
| |
| // this was the last point we needed, we can compute the derivatives |
| setStepSize(t[t.length - 1].subtract(t[0]).divide(t.length - 1)); |
| |
| // first scaled derivative |
| scaled = MathArrays.buildArray(getField(), yDot[0].length); |
| for (int j = 0; j < scaled.length; ++j) { |
| scaled[j] = yDot[0][j].multiply(getStepSize()); |
| } |
| |
| // higher order derivatives |
| nordsieck = initializeHighOrderDerivatives(getStepSize(), t, y, yDot); |
| |
| // stop the integrator now that all needed steps have been handled |
| setStepStart(savedStart); |
| throw new InitializationCompletedMarkerException(); |
| |
| } |
| |
| } |
| |
| /** {@inheritDoc} */ |
| @Override |
| public void init(final FieldODEStateAndDerivative<T> initialState, T finalTime) { |
| // nothing to do |
| } |
| |
| } |
| |
| /** Marker exception used ONLY to stop the starter integrator after first step. */ |
| private static class InitializationCompletedMarkerException |
| extends RuntimeException { |
| |
| /** Serializable version identifier. */ |
| private static final long serialVersionUID = -1914085471038046418L; |
| |
| /** Simple constructor. */ |
| InitializationCompletedMarkerException() { |
| super((Throwable) null); |
| } |
| |
| } |
| |
| } |