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.math3.ode.nonstiff;
019
020import org.apache.commons.math3.Field;
021import org.apache.commons.math3.RealFieldElement;
022import org.apache.commons.math3.exception.DimensionMismatchException;
023import org.apache.commons.math3.exception.MaxCountExceededException;
024import org.apache.commons.math3.exception.NumberIsTooSmallException;
025import org.apache.commons.math3.exception.util.LocalizedFormats;
026import org.apache.commons.math3.ode.AbstractFieldIntegrator;
027import org.apache.commons.math3.ode.FieldEquationsMapper;
028import org.apache.commons.math3.ode.FieldODEState;
029import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
030import org.apache.commons.math3.util.FastMath;
031import org.apache.commons.math3.util.MathArrays;
032import org.apache.commons.math3.util.MathUtils;
033
034/**
035 * This abstract class holds the common part of all adaptive
036 * stepsize integrators for Ordinary Differential Equations.
037 *
038 * <p>These algorithms perform integration with stepsize control, which
039 * means the user does not specify the integration step but rather a
040 * tolerance on error. The error threshold is computed as
041 * <pre>
042 * threshold_i = absTol_i + relTol_i * max (abs (ym), abs (ym+1))
043 * </pre>
044 * where absTol_i is the absolute tolerance for component i of the
045 * state vector and relTol_i is the relative tolerance for the same
046 * component. The user can also use only two scalar values absTol and
047 * relTol which will be used for all components.
048 * </p>
049 * <p>
050 * Note that <em>only</em> the {@link FieldODEState#getState() main part}
051 * of the state vector is used for stepsize control. The {@link
052 * FieldODEState#getSecondaryState(int) secondary parts} of the state
053 * vector are explicitly ignored for stepsize control.
054 * </p>
055 *
056 * <p>If the estimated error for ym+1 is such that
057 * <pre>
058 * sqrt((sum (errEst_i / threshold_i)^2 ) / n) < 1
059 * </pre>
060 *
061 * (where n is the main set dimension) then the step is accepted,
062 * otherwise the step is rejected and a new attempt is made with a new
063 * stepsize.</p>
064 *
065 * @param <T> the type of the field elements
066 * @since 3.6
067 *
068 */
069
070public abstract class AdaptiveStepsizeFieldIntegrator<T extends RealFieldElement<T>>
071    extends AbstractFieldIntegrator<T> {
072
073    /** Allowed absolute scalar error. */
074    protected double scalAbsoluteTolerance;
075
076    /** Allowed relative scalar error. */
077    protected double scalRelativeTolerance;
078
079    /** Allowed absolute vectorial error. */
080    protected double[] vecAbsoluteTolerance;
081
082    /** Allowed relative vectorial error. */
083    protected double[] vecRelativeTolerance;
084
085    /** Main set dimension. */
086    protected int mainSetDimension;
087
088    /** User supplied initial step. */
089    private T initialStep;
090
091    /** Minimal step. */
092    private T minStep;
093
094    /** Maximal step. */
095    private T maxStep;
096
097    /** Build an integrator with the given stepsize bounds.
098     * The default step handler does nothing.
099     * @param field field to which the time and state vector elements belong
100     * @param name name of the method
101     * @param minStep minimal step (sign is irrelevant, regardless of
102     * integration direction, forward or backward), the last step can
103     * be smaller than this
104     * @param maxStep maximal step (sign is irrelevant, regardless of
105     * integration direction, forward or backward), the last step can
106     * be smaller than this
107     * @param scalAbsoluteTolerance allowed absolute error
108     * @param scalRelativeTolerance allowed relative error
109     */
110    public AdaptiveStepsizeFieldIntegrator(final Field<T> field, final String name,
111                                           final double minStep, final double maxStep,
112                                           final double scalAbsoluteTolerance,
113                                           final double scalRelativeTolerance) {
114
115        super(field, name);
116        setStepSizeControl(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
117        resetInternalState();
118
119    }
120
121    /** Build an integrator with the given stepsize bounds.
122     * The default step handler does nothing.
123     * @param field field to which the time and state vector elements belong
124     * @param name name of the method
125     * @param minStep minimal step (sign is irrelevant, regardless of
126     * integration direction, forward or backward), the last step can
127     * be smaller than this
128     * @param maxStep maximal step (sign is irrelevant, regardless of
129     * integration direction, forward or backward), the last step can
130     * be smaller than this
131     * @param vecAbsoluteTolerance allowed absolute error
132     * @param vecRelativeTolerance allowed relative error
133     */
134    public AdaptiveStepsizeFieldIntegrator(final Field<T> field, final String name,
135                                           final double minStep, final double maxStep,
136                                           final double[] vecAbsoluteTolerance,
137                                           final double[] vecRelativeTolerance) {
138
139        super(field, name);
140        setStepSizeControl(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
141        resetInternalState();
142
143    }
144
145    /** Set the adaptive step size control parameters.
146     * <p>
147     * A side effect of this method is to also reset the initial
148     * step so it will be automatically computed by the integrator
149     * if {@link #setInitialStepSize(RealFieldElement) setInitialStepSize}
150     * is not called by the user.
151     * </p>
152     * @param minimalStep minimal step (must be positive even for backward
153     * integration), the last step can be smaller than this
154     * @param maximalStep maximal step (must be positive even for backward
155     * integration)
156     * @param absoluteTolerance allowed absolute error
157     * @param relativeTolerance allowed relative error
158     */
159    public void setStepSizeControl(final double minimalStep, final double maximalStep,
160                                   final double absoluteTolerance,
161                                   final double relativeTolerance) {
162
163        minStep     = getField().getZero().add(FastMath.abs(minimalStep));
164        maxStep     = getField().getZero().add(FastMath.abs(maximalStep));
165        initialStep = getField().getOne().negate();
166
167        scalAbsoluteTolerance = absoluteTolerance;
168        scalRelativeTolerance = relativeTolerance;
169        vecAbsoluteTolerance  = null;
170        vecRelativeTolerance  = null;
171
172    }
173
174    /** Set the adaptive step size control parameters.
175     * <p>
176     * A side effect of this method is to also reset the initial
177     * step so it will be automatically computed by the integrator
178     * if {@link #setInitialStepSize(RealFieldElement) setInitialStepSize}
179     * is not called by the user.
180     * </p>
181     * @param minimalStep minimal step (must be positive even for backward
182     * integration), the last step can be smaller than this
183     * @param maximalStep maximal step (must be positive even for backward
184     * integration)
185     * @param absoluteTolerance allowed absolute error
186     * @param relativeTolerance allowed relative error
187     */
188    public void setStepSizeControl(final double minimalStep, final double maximalStep,
189                                   final double[] absoluteTolerance,
190                                   final double[] relativeTolerance) {
191
192        minStep     = getField().getZero().add(FastMath.abs(minimalStep));
193        maxStep     = getField().getZero().add(FastMath.abs(maximalStep));
194        initialStep = getField().getOne().negate();
195
196        scalAbsoluteTolerance = 0;
197        scalRelativeTolerance = 0;
198        vecAbsoluteTolerance  = absoluteTolerance.clone();
199        vecRelativeTolerance  = relativeTolerance.clone();
200
201    }
202
203    /** Set the initial step size.
204     * <p>This method allows the user to specify an initial positive
205     * step size instead of letting the integrator guess it by
206     * itself. If this method is not called before integration is
207     * started, the initial step size will be estimated by the
208     * integrator.</p>
209     * @param initialStepSize initial step size to use (must be positive even
210     * for backward integration ; providing a negative value or a value
211     * outside of the min/max step interval will lead the integrator to
212     * ignore the value and compute the initial step size by itself)
213     */
214    public void setInitialStepSize(final T initialStepSize) {
215        if (initialStepSize.subtract(minStep).getReal() < 0 ||
216            initialStepSize.subtract(maxStep).getReal() > 0) {
217            initialStep = getField().getOne().negate();
218        } else {
219            initialStep = initialStepSize;
220        }
221    }
222
223    /** {@inheritDoc} */
224    @Override
225    protected void sanityChecks(final FieldODEState<T> eqn, final T t)
226        throws DimensionMismatchException, NumberIsTooSmallException {
227
228        super.sanityChecks(eqn, t);
229
230        mainSetDimension = eqn.getStateDimension();
231
232        if (vecAbsoluteTolerance != null && vecAbsoluteTolerance.length != mainSetDimension) {
233            throw new DimensionMismatchException(mainSetDimension, vecAbsoluteTolerance.length);
234        }
235
236        if (vecRelativeTolerance != null && vecRelativeTolerance.length != mainSetDimension) {
237            throw new DimensionMismatchException(mainSetDimension, vecRelativeTolerance.length);
238        }
239
240    }
241
242    /** Initialize the integration step.
243     * @param forward forward integration indicator
244     * @param order order of the method
245     * @param scale scaling vector for the state vector (can be shorter than state vector)
246     * @param state0 state at integration start time
247     * @param mapper mapper for all the equations
248     * @return first integration step
249     * @exception MaxCountExceededException if the number of functions evaluations is exceeded
250     * @exception DimensionMismatchException if arrays dimensions do not match equations settings
251     */
252    public T initializeStep(final boolean forward, final int order, final T[] scale,
253                            final FieldODEStateAndDerivative<T> state0,
254                            final FieldEquationsMapper<T> mapper)
255        throws MaxCountExceededException, DimensionMismatchException {
256
257        if (initialStep.getReal() > 0) {
258            // use the user provided value
259            return forward ? initialStep : initialStep.negate();
260        }
261
262        // very rough first guess : h = 0.01 * ||y/scale|| / ||y'/scale||
263        // this guess will be used to perform an Euler step
264        final T[] y0    = mapper.mapState(state0);
265        final T[] yDot0 = mapper.mapDerivative(state0);
266        T yOnScale2    = getField().getZero();
267        T yDotOnScale2 = getField().getZero();
268        for (int j = 0; j < scale.length; ++j) {
269            final T ratio    = y0[j].divide(scale[j]);
270            yOnScale2        = yOnScale2.add(ratio.multiply(ratio));
271            final T ratioDot = yDot0[j].divide(scale[j]);
272            yDotOnScale2     = yDotOnScale2.add(ratioDot.multiply(ratioDot));
273        }
274
275        T h = (yOnScale2.getReal() < 1.0e-10 || yDotOnScale2.getReal() < 1.0e-10) ?
276              getField().getZero().add(1.0e-6) :
277              yOnScale2.divide(yDotOnScale2).sqrt().multiply(0.01);
278        if (! forward) {
279            h = h.negate();
280        }
281
282        // perform an Euler step using the preceding rough guess
283        final T[] y1 = MathArrays.buildArray(getField(), y0.length);
284        for (int j = 0; j < y0.length; ++j) {
285            y1[j] = y0[j].add(yDot0[j].multiply(h));
286        }
287        final T[] yDot1 = computeDerivatives(state0.getTime().add(h), y1);
288
289        // estimate the second derivative of the solution
290        T yDDotOnScale = getField().getZero();
291        for (int j = 0; j < scale.length; ++j) {
292            final T ratioDotDot = yDot1[j].subtract(yDot0[j]).divide(scale[j]);
293            yDDotOnScale = yDDotOnScale.add(ratioDotDot.multiply(ratioDotDot));
294        }
295        yDDotOnScale = yDDotOnScale.sqrt().divide(h);
296
297        // step size is computed such that
298        // h^order * max (||y'/tol||, ||y''/tol||) = 0.01
299        final T maxInv2 = MathUtils.max(yDotOnScale2.sqrt(), yDDotOnScale);
300        final T h1 = maxInv2.getReal() < 1.0e-15 ?
301                     MathUtils.max(getField().getZero().add(1.0e-6), h.abs().multiply(0.001)) :
302                     maxInv2.multiply(100).reciprocal().pow(1.0 / order);
303        h = MathUtils.min(h.abs().multiply(100), h1);
304        h = MathUtils.max(h, state0.getTime().abs().multiply(1.0e-12));  // avoids cancellation when computing t1 - t0
305        h = MathUtils.max(minStep, MathUtils.min(maxStep, h));
306        if (! forward) {
307            h = h.negate();
308        }
309
310        return h;
311
312    }
313
314    /** Filter the integration step.
315     * @param h signed step
316     * @param forward forward integration indicator
317     * @param acceptSmall if true, steps smaller than the minimal value
318     * are silently increased up to this value, if false such small
319     * steps generate an exception
320     * @return a bounded integration step (h if no bound is reach, or a bounded value)
321     * @exception NumberIsTooSmallException if the step is too small and acceptSmall is false
322     */
323    protected T filterStep(final T h, final boolean forward, final boolean acceptSmall)
324        throws NumberIsTooSmallException {
325
326        T filteredH = h;
327        if (h.abs().subtract(minStep).getReal() < 0) {
328            if (acceptSmall) {
329                filteredH = forward ? minStep : minStep.negate();
330            } else {
331                throw new NumberIsTooSmallException(LocalizedFormats.MINIMAL_STEPSIZE_REACHED_DURING_INTEGRATION,
332                                                    h.abs().getReal(), minStep.getReal(), true);
333            }
334        }
335
336        if (filteredH.subtract(maxStep).getReal() > 0) {
337            filteredH = maxStep;
338        } else if (filteredH.add(maxStep).getReal() < 0) {
339            filteredH = maxStep.negate();
340        }
341
342        return filteredH;
343
344    }
345
346    /** Reset internal state to dummy values. */
347    protected void resetInternalState() {
348        setStepStart(null);
349        setStepSize(minStep.multiply(maxStep).sqrt());
350    }
351
352    /** Get the minimal step.
353     * @return minimal step
354     */
355    public T getMinStep() {
356        return minStep;
357    }
358
359    /** Get the maximal step.
360     * @return maximal step
361     */
362    public T getMaxStep() {
363        return maxStep;
364    }
365
366}