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