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.CalculusFieldElement;
21  import org.hipparchus.exception.LocalizedCoreFormats;
22  import org.hipparchus.exception.MathIllegalArgumentException;
23  import org.hipparchus.ode.LocalizedODEFormats;
24  import org.hipparchus.util.FastMath;
25  
26  /** Helper for adaptive stepsize control.
27   * @since 2.0
28   */
29  
30  public class StepsizeHelper {
31  
32      /** Allowed absolute scalar error. */
33      private double scalAbsoluteTolerance;
34  
35      /** Allowed relative scalar error. */
36      private double scalRelativeTolerance;
37  
38      /** Allowed absolute vectorial error. */
39      private double[] vecAbsoluteTolerance;
40  
41      /** Allowed relative vectorial error. */
42      private double[] vecRelativeTolerance;
43  
44      /** Main set dimension. */
45      private int mainSetDimension;
46  
47      /** User supplied initial step. */
48      private double initialStep;
49  
50      /** Minimal step. */
51      private double minStep;
52  
53      /** Maximal step. */
54      private double maxStep;
55  
56      /** Simple constructor.
57       * @param minStep minimal step (sign is irrelevant, regardless of
58       * integration direction, forward or backward), the last step can
59       * be smaller than this
60       * @param maxStep maximal step (sign is irrelevant, regardless of
61       * integration direction, forward or backward), the last step can
62       * be smaller than this
63       * @param scalAbsoluteTolerance allowed absolute error
64       * @param scalRelativeTolerance allowed relative error
65       */
66      public StepsizeHelper(final double minStep, final double maxStep,
67                            final double scalAbsoluteTolerance,
68                            final double scalRelativeTolerance) {
69          this.minStep     = FastMath.abs(minStep);
70          this.maxStep     = FastMath.abs(maxStep);
71          this.initialStep = -1;
72  
73          this.scalAbsoluteTolerance = scalAbsoluteTolerance;
74          this.scalRelativeTolerance = scalRelativeTolerance;
75          this.vecAbsoluteTolerance  = null;
76          this.vecRelativeTolerance  = null;
77      }
78  
79      /** Simple constructor..
80       * @param minStep minimal step (sign is irrelevant, regardless of
81       * integration direction, forward or backward), the last step can
82       * be smaller than this
83       * @param maxStep maximal step (sign is irrelevant, regardless of
84       * integration direction, forward or backward), the last step can
85       * be smaller than this
86       * @param vecAbsoluteTolerance allowed absolute error
87       * @param vecRelativeTolerance allowed relative error
88       */
89      public StepsizeHelper(final double minStep, final double maxStep,
90                            final double[] vecAbsoluteTolerance,
91                            final double[] vecRelativeTolerance) {
92  
93          this.minStep     = FastMath.abs(minStep);
94          this.maxStep     = FastMath.abs(maxStep);
95          this.initialStep = -1;
96  
97         this.scalAbsoluteTolerance = 0;
98         this.scalRelativeTolerance = 0;
99         this.vecAbsoluteTolerance  = vecAbsoluteTolerance.clone();
100        this.vecRelativeTolerance  = vecRelativeTolerance.clone();
101 
102     }
103 
104     /** Set main set dimension.
105      * @param mainSetDimension dimension of the main set
106      * @exception MathIllegalArgumentException if adaptive step size integrators
107      * tolerance arrays dimensions are not compatible with equations settings
108      */
109     protected void setMainSetDimension(final int mainSetDimension) throws MathIllegalArgumentException {
110         this.mainSetDimension = mainSetDimension;
111 
112         if (vecAbsoluteTolerance != null && vecAbsoluteTolerance.length != mainSetDimension) {
113             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
114                                                    mainSetDimension, vecAbsoluteTolerance.length);
115         }
116 
117         if (vecRelativeTolerance != null && vecRelativeTolerance.length != mainSetDimension) {
118             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
119                                                    mainSetDimension, vecRelativeTolerance.length);
120         }
121     }
122 
123     /** Get the main set dimension.
124      * @return main set dimension
125      */
126     public int getMainSetDimension() {
127         return mainSetDimension;
128     }
129 
130     /** Get the relative tolerance for one component.
131      * @param i component to select
132      * @return relative tolerance for selected component
133      */
134     public double getRelativeTolerance(final int i) {
135         return vecAbsoluteTolerance == null ? scalRelativeTolerance : vecRelativeTolerance[i];
136     }
137 
138     /** Get the tolerance for one component.
139      * @param i component to select
140      * @param scale scale factor for relative tolerance (i.e. y[i])
141      * @return tolerance for selected component
142      */
143     public double getTolerance(final int i, final double scale) {
144         return vecAbsoluteTolerance == null ?
145                scalAbsoluteTolerance   + scalRelativeTolerance   * scale :
146                vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * scale;
147     }
148 
149     /** Get the tolerance for one component.
150      * @param i component to select
151      * @param scale scale factor for relative tolerance (i.e. y[i])
152      * @param <T> type of the field elements
153      * @return tolerance for selected component
154      */
155     public <T extends CalculusFieldElement<T>> T getTolerance(final int i, final T scale) {
156         return vecAbsoluteTolerance == null ?
157                scale.multiply(scalRelativeTolerance).add(scalAbsoluteTolerance) :
158                scale.multiply(vecRelativeTolerance[i]).add(vecAbsoluteTolerance[i]);
159     }
160 
161     /** Filter the integration step.
162      * @param h signed step
163      * @param forward forward integration indicator
164      * @param acceptSmall if true, steps smaller than the minimal value
165      * are silently increased up to this value, if false such small
166      * steps generate an exception
167      * @return a bounded integration step (h if no bound is reach, or a bounded value)
168      * @exception MathIllegalArgumentException if the step is too small and acceptSmall is false
169      */
170     public double filterStep(final double h, final boolean forward, final boolean acceptSmall)
171         throws MathIllegalArgumentException {
172 
173         double filteredH = h;
174         if (FastMath.abs(h) < minStep) {
175             if (acceptSmall) {
176                 filteredH = forward ? minStep : -minStep;
177             } else {
178                 throw new MathIllegalArgumentException(LocalizedODEFormats.MINIMAL_STEPSIZE_REACHED_DURING_INTEGRATION,
179                                                        FastMath.abs(h), minStep, true);
180             }
181         }
182 
183         if (filteredH > maxStep) {
184             filteredH = maxStep;
185         } else if (filteredH < -maxStep) {
186             filteredH = -maxStep;
187         }
188 
189         return filteredH;
190 
191     }
192 
193     /** Filter the integration step.
194      * @param h signed step
195      * @param forward forward integration indicator
196      * @param acceptSmall if true, steps smaller than the minimal value
197      * are silently increased up to this value, if false such small
198      * steps generate an exception
199      * @param <T> type of the field elements
200      * @return a bounded integration step (h if no bound is reach, or a bounded value)
201      * @exception MathIllegalArgumentException if the step is too small and acceptSmall is false
202      */
203     public <T extends CalculusFieldElement<T>> T filterStep(final T h, final boolean forward, final boolean acceptSmall)
204         throws MathIllegalArgumentException {
205 
206         T filteredH = h;
207         if (h.abs().subtract(minStep).getReal() < 0) {
208             if (acceptSmall) {
209                 filteredH = h.getField().getZero().add(forward ? minStep : -minStep);
210             } else {
211                 throw new MathIllegalArgumentException(LocalizedODEFormats.MINIMAL_STEPSIZE_REACHED_DURING_INTEGRATION,
212                                                        FastMath.abs(h.getReal()), minStep, true);
213             }
214         }
215 
216         if (filteredH.subtract(maxStep).getReal() > 0) {
217             filteredH = h.getField().getZero().add(maxStep);
218         } else if (filteredH.add(maxStep).getReal() < 0) {
219             filteredH = h.getField().getZero().add(-maxStep);
220         }
221 
222         return filteredH;
223 
224     }
225 
226     /** Set the initial step size.
227      * <p>This method allows the user to specify an initial positive
228      * step size instead of letting the integrator guess it by
229      * itself. If this method is not called before integration is
230      * started, the initial step size will be estimated by the
231      * integrator.</p>
232      * @param initialStepSize initial step size to use (must be positive even
233      * for backward integration ; providing a negative value or a value
234      * outside of the min/max step interval will lead the integrator to
235      * ignore the value and compute the initial step size by itself)
236      */
237     public void setInitialStepSize(final double initialStepSize) {
238         if ((initialStepSize < minStep) || (initialStepSize > maxStep)) {
239             initialStep = -1.0;
240         } else {
241             initialStep = initialStepSize;
242         }
243     }
244 
245     /** Get the initial step.
246      * @return initial step
247      */
248     public double getInitialStep() {
249         return initialStep;
250     }
251 
252     /** Get the minimal step.
253      * @return minimal step
254      */
255     public double getMinStep() {
256         return minStep;
257     }
258 
259     /** Get the maximal step.
260      * @return maximal step
261      */
262     public double getMaxStep() {
263         return maxStep;
264     }
265 
266     /** Get a dummy step size.
267      * @return geometric mean of {@link #getMinStep()} and {@link #getMaxStep()}
268      */
269     public double getDummyStepsize() {
270         return FastMath.sqrt(minStep * maxStep);
271     }
272 
273 }