View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  
23  package org.hipparchus.ode.nonstiff;
24  
25  import org.hipparchus.CalculusFieldElement;
26  import org.hipparchus.Field;
27  import org.hipparchus.exception.MathIllegalArgumentException;
28  import org.hipparchus.exception.MathIllegalStateException;
29  import org.hipparchus.ode.AbstractFieldIntegrator;
30  import org.hipparchus.ode.FieldEquationsMapper;
31  import org.hipparchus.ode.FieldODEState;
32  import org.hipparchus.ode.FieldODEStateAndDerivative;
33  import org.hipparchus.util.FastMath;
34  import org.hipparchus.util.MathArrays;
35  
36  /**
37   * This abstract class holds the common part of all adaptive
38   * stepsize integrators for Ordinary Differential Equations.
39   *
40   * <p>These algorithms perform integration with stepsize control, which
41   * means the user does not specify the integration step but rather a
42   * tolerance on error. The error threshold is computed as
43   * </p>
44   * <pre>
45   * threshold_i = absTol_i + relTol_i * max (abs (ym), abs (ym+1))
46   * </pre>
47   * <p>
48   * where absTol_i is the absolute tolerance for component i of the
49   * state vector and relTol_i is the relative tolerance for the same
50   * component. The user can also use only two scalar values absTol and
51   * relTol which will be used for all components.
52   * </p>
53   * <p>
54   * Note that <em>only</em> the {@link FieldODEState#getPrimaryState() main part}
55   * of the state vector is used for stepsize control. The {@link
56   * FieldODEState#getSecondaryState(int) secondary parts} of the state
57   * vector are explicitly ignored for stepsize control.
58   * </p>
59   *
60   * <p>If the estimated error for ym+1 is such that</p>
61   * <pre>
62   * sqrt((sum (errEst_i / threshold_i)^2 ) / n) &lt; 1
63   * </pre>
64   *
65   * <p>(where n is the main set dimension) then the step is accepted,
66   * otherwise the step is rejected and a new attempt is made with a new
67   * stepsize.</p>
68   *
69   * @param <T> the type of the field elements
70   *
71   */
72  
73  public abstract class AdaptiveStepsizeFieldIntegrator<T extends CalculusFieldElement<T>>
74      extends AbstractFieldIntegrator<T> {
75  
76      /** Helper for step size control. */
77      private StepsizeHelper stepsizeHelper;
78  
79      /** Build an integrator with the given stepsize bounds.
80       * The default step handler does nothing.
81       * @param field field to which the time and state vector elements belong
82       * @param name name of the method
83       * @param minStep minimal step (sign is irrelevant, regardless of
84       * integration direction, forward or backward), the last step can
85       * be smaller than this
86       * @param maxStep maximal step (sign is irrelevant, regardless of
87       * integration direction, forward or backward), the last step can
88       * be smaller than this
89       * @param scalAbsoluteTolerance allowed absolute error
90       * @param scalRelativeTolerance allowed relative error
91       */
92      public AdaptiveStepsizeFieldIntegrator(final Field<T> field, final String name,
93                                             final double minStep, final double maxStep,
94                                             final double scalAbsoluteTolerance,
95                                             final double scalRelativeTolerance) {
96          super(field, name);
97          stepsizeHelper = new StepsizeHelper(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
98          resetInternalState();
99      }
100 
101     /** Build an integrator with the given stepsize bounds.
102      * The default step handler does nothing.
103      * @param field field to which the time and state vector elements belong
104      * @param name name of the method
105      * @param minStep minimal step (sign is irrelevant, regardless of
106      * integration direction, forward or backward), the last step can
107      * be smaller than this
108      * @param maxStep maximal step (sign is irrelevant, regardless of
109      * integration direction, forward or backward), the last step can
110      * be smaller than this
111      * @param vecAbsoluteTolerance allowed absolute error
112      * @param vecRelativeTolerance allowed relative error
113      */
114     public AdaptiveStepsizeFieldIntegrator(final Field<T> field, final String name,
115                                            final double minStep, final double maxStep,
116                                            final double[] vecAbsoluteTolerance,
117                                            final double[] vecRelativeTolerance) {
118         super(field, name);
119         stepsizeHelper = new StepsizeHelper(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
120         resetInternalState();
121     }
122 
123     /** Set the adaptive step size control parameters.
124      * <p>
125      * A side effect of this method is to also reset the initial
126      * step so it will be automatically computed by the integrator
127      * if {@link #setInitialStepSize(double) setInitialStepSize}
128      * is not called by the user.
129      * </p>
130      * @param minimalStep minimal step (must be positive even for backward
131      * integration), the last step can be smaller than this
132      * @param maximalStep maximal step (must be positive even for backward
133      * integration)
134      * @param absoluteTolerance allowed absolute error
135      * @param relativeTolerance allowed relative error
136      */
137     public void setStepSizeControl(final double minimalStep, final double maximalStep,
138                                    final double absoluteTolerance, final double relativeTolerance) {
139         stepsizeHelper = new StepsizeHelper(minimalStep, maximalStep, absoluteTolerance, relativeTolerance);
140     }
141 
142     /** Set the adaptive step size control parameters.
143      * <p>
144      * A side effect of this method is to also reset the initial
145      * step so it will be automatically computed by the integrator
146      * if {@link #setInitialStepSize(double) setInitialStepSize}
147      * is not called by the user.
148      * </p>
149      * @param minimalStep minimal step (must be positive even for backward
150      * integration), the last step can be smaller than this
151      * @param maximalStep maximal step (must be positive even for backward
152      * integration)
153      * @param absoluteTolerance allowed absolute error
154      * @param relativeTolerance allowed relative error
155      */
156     public void setStepSizeControl(final double minimalStep, final double maximalStep,
157                                    final double[] absoluteTolerance,
158                                    final double[] relativeTolerance) {
159         stepsizeHelper = new StepsizeHelper(minimalStep, maximalStep, absoluteTolerance, relativeTolerance);
160     }
161 
162     /** Get the stepsize helper.
163      * @return stepsize helper
164      * @since 2.0
165      */
166     protected StepsizeHelper getStepSizeHelper() {
167         return stepsizeHelper;
168     }
169 
170     /** Set the initial step size.
171      * <p>This method allows the user to specify an initial positive
172      * step size instead of letting the integrator guess it by
173      * itself. If this method is not called before integration is
174      * started, the initial step size will be estimated by the
175      * integrator.</p>
176      * @param initialStepSize initial step size to use (must be positive even
177      * for backward integration ; providing a negative value or a value
178      * outside of the min/max step interval will lead the integrator to
179      * ignore the value and compute the initial step size by itself)
180      */
181     public void setInitialStepSize(final double initialStepSize) {
182         stepsizeHelper.setInitialStepSize(initialStepSize);
183     }
184 
185     /** {@inheritDoc} */
186     @Override
187     protected void sanityChecks(final FieldODEState<T> initialState, final T t)
188         throws MathIllegalArgumentException {
189         super.sanityChecks(initialState, t);
190         stepsizeHelper.setMainSetDimension(initialState.getPrimaryStateDimension());
191     }
192 
193     /** Initialize the integration step.
194      * @param forward forward integration indicator
195      * @param order order of the method
196      * @param scale scaling vector for the state vector (can be shorter than state vector)
197      * @param state0 state at integration start time
198      * @param mapper mapper for all the equations
199      * @return first integration step
200      * @exception MathIllegalStateException if the number of functions evaluations is exceeded
201      * @exception MathIllegalArgumentException if arrays dimensions do not match equations settings
202      */
203     public double initializeStep(final boolean forward, final int order, final T[] scale,
204                                  final FieldODEStateAndDerivative<T> state0,
205                                  final FieldEquationsMapper<T> mapper)
206         throws MathIllegalArgumentException, MathIllegalStateException {
207 
208         if (stepsizeHelper.getInitialStep() > 0) {
209             // use the user provided value
210             return forward ? stepsizeHelper.getInitialStep() : -stepsizeHelper.getInitialStep();
211         }
212 
213         // very rough first guess : h = 0.01 * ||y/scale|| / ||y'/scale||
214         // this guess will be used to perform an Euler step
215         final T[] y0    = state0.getCompleteState();
216         final T[] yDot0 = state0.getCompleteDerivative();
217         double yOnScale2     = 0;
218         double yDotOnScale2  = 0;
219         for (int j = 0; j < scale.length; ++j) {
220             final double ratio    = y0[j].getReal() / scale[j].getReal();
221             yOnScale2            += ratio * ratio;
222             final double ratioDot = yDot0[j].getReal() / scale[j].getReal();
223             yDotOnScale2         += ratioDot * ratioDot;
224         }
225 
226         double h = ((yOnScale2 < 1.0e-10) || (yDotOnScale2 < 1.0e-10)) ?
227                    1.0e-6 : (0.01 * FastMath.sqrt(yOnScale2 / yDotOnScale2));
228         if (h > getMaxStep()) {
229             h = getMaxStep();
230         }
231         if (! forward) {
232             h = -h;
233         }
234 
235         // perform an Euler step using the preceding rough guess
236         final T[] y1 = MathArrays.buildArray(getField(), y0.length);
237         for (int j = 0; j < y0.length; ++j) {
238             y1[j] = y0[j].add(yDot0[j].multiply(h));
239         }
240         final T[] yDot1 = computeDerivatives(state0.getTime().add(h), y1);
241 
242         // estimate the second derivative of the solution
243         double yDDotOnScale = 0;
244         for (int j = 0; j < scale.length; ++j) {
245             final double ratioDotDot = (yDot1[j].getReal() - yDot0[j].getReal()) / scale[j].getReal();
246             yDDotOnScale += ratioDotDot * ratioDotDot;
247         }
248         yDDotOnScale = FastMath.sqrt(yDDotOnScale) / h;
249 
250         // step size is computed such that
251         // h^order * max (||y'/tol||, ||y''/tol||) = 0.01
252         final double maxInv2 = FastMath.max(FastMath.sqrt(yDotOnScale2), yDDotOnScale);
253         final double h1 = (maxInv2 < 1.0e-15) ?
254                           FastMath.max(1.0e-6, 0.001 * FastMath.abs(h)) :
255                           FastMath.pow(0.01 / maxInv2, 1.0 / order);
256         h = FastMath.min(100.0 * FastMath.abs(h), h1);
257         h = FastMath.max(h, 1.0e-12 * FastMath.abs(state0.getTime().getReal()));  // avoids cancellation when computing t1 - t0
258         if (h < getMinStep()) {
259             h = getMinStep();
260         }
261         if (h > getMaxStep()) {
262             h = getMaxStep();
263         }
264 
265         if (! forward) {
266             h = -h;
267         }
268 
269         return h;
270 
271     }
272 
273     /** Reset internal state to dummy values. */
274     protected void resetInternalState() {
275         setStepStart(null);
276         setStepSize(getField().getZero().add(stepsizeHelper.getDummyStepsize()));
277     }
278 
279     /** Get the minimal step.
280      * @return minimal step
281      */
282     public double getMinStep() {
283         return stepsizeHelper.getMinStep();
284     }
285 
286     /** Get the maximal step.
287      * @return maximal step
288      */
289     public double getMaxStep() {
290         return stepsizeHelper.getMaxStep();
291     }
292 
293 }