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().newInstance(maxStep);
218 } else if (filteredH.add(maxStep).getReal() < 0) {
219 filteredH = h.getField().getZero().newInstance(-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 }