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.FieldEquationsMapper;
30 import org.hipparchus.ode.FieldExpandableODE;
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 import org.hipparchus.util.MathUtils;
36
37 /**
38 * This class implements the common part of all embedded Runge-Kutta
39 * integrators for Ordinary Differential Equations.
40 *
41 * <p>These methods are embedded explicit Runge-Kutta methods with two
42 * sets of coefficients allowing to estimate the error, their Butcher
43 * arrays are as follows :</p>
44 * <pre>
45 * 0 |
46 * c2 | a21
47 * c3 | a31 a32
48 * ... | ...
49 * cs | as1 as2 ... ass-1
50 * |--------------------------
51 * | b1 b2 ... bs-1 bs
52 * | b'1 b'2 ... b's-1 b's
53 * </pre>
54 *
55 * <p>In fact, we rather use the array defined by ej = bj - b'j to
56 * compute directly the error rather than computing two estimates and
57 * then comparing them.</p>
58 *
59 * <p>Some methods are qualified as <i>fsal</i> (first same as last)
60 * methods. This means the last evaluation of the derivatives in one
61 * step is the same as the first in the next step. Then, this
62 * evaluation can be reused from one step to the next one and the cost
63 * of such a method is really s-1 evaluations despite the method still
64 * has s stages. This behaviour is true only for successful steps, if
65 * the step is rejected after the error estimation phase, no
66 * evaluation is saved. For an <i>fsal</i> method, we have cs = 1 and
67 * asi = bi for all i.</p>
68 *
69 * @param <T> the type of the field elements
70 */
71
72 public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends CalculusFieldElement<T>>
73 extends AdaptiveStepsizeFieldIntegrator<T>
74 implements FieldExplicitRungeKuttaIntegrator<T> {
75
76 /** Index of the pre-computed derivative for <i>fsal</i> methods. */
77 private final int fsal;
78
79 /** Time steps from Butcher array (without the first zero). */
80 private final T[] c;
81
82 /** Internal weights from Butcher array (without the first empty row). */
83 private final T[][] a;
84
85 /** External weights for the high order method from Butcher array. */
86 private final T[] b;
87
88 /** Time steps from Butcher array (without the first zero). */
89 private double[] realC = new double[0];
90
91 /** Internal weights from Butcher array (without the first empty row). Real version, optional. */
92 private double[][] realA = new double[0][];
93
94 /** External weights for the high order method from Butcher array. Real version, optional. */
95 private double[] realB = new double[0];
96
97 /** Stepsize control exponent. */
98 private final double exp;
99
100 /** Safety factor for stepsize control. */
101 private T safety;
102
103 /** Minimal reduction factor for stepsize control. */
104 private T minReduction;
105
106 /** Maximal growth factor for stepsize control. */
107 private T maxGrowth;
108
109 /** Flag setting whether coefficients in Butcher array are interpreted as Field or real numbers. */
110 private boolean usingFieldCoefficients;
111
112 /** Build a Runge-Kutta integrator with the given Butcher array.
113 * @param field field to which the time and state vector elements belong
114 * @param name name of the method
115 * @param fsal index of the pre-computed derivative for <i>fsal</i> methods
116 * or -1 if method is not <i>fsal</i>
117 * @param minStep minimal step (sign is irrelevant, regardless of
118 * integration direction, forward or backward), the last step can
119 * be smaller than this
120 * @param maxStep maximal step (sign is irrelevant, regardless of
121 * integration direction, forward or backward), the last step can
122 * be smaller than this
123 * @param scalAbsoluteTolerance allowed absolute error
124 * @param scalRelativeTolerance allowed relative error
125 */
126 protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal,
127 final double minStep, final double maxStep,
128 final double scalAbsoluteTolerance,
129 final double scalRelativeTolerance) {
130
131 super(field, name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
132
133 this.fsal = fsal;
134 this.c = getC();
135 this.a = getA();
136 this.b = getB();
137
138 exp = -1.0 / getOrder();
139
140 // set the default values of the algorithm control parameters
141 setSafety(field.getZero().add(0.9));
142 setMinReduction(field.getZero().add(0.2));
143 setMaxGrowth(field.getZero().add(10.0));
144
145 }
146
147 /** Build a Runge-Kutta integrator with the given Butcher array.
148 * @param field field to which the time and state vector elements belong
149 * @param name name of the method
150 * @param fsal index of the pre-computed derivative for <i>fsal</i> methods
151 * or -1 if method is not <i>fsal</i>
152 * @param minStep minimal step (must be positive even for backward
153 * integration), the last step can be smaller than this
154 * @param maxStep maximal step (must be positive even for backward
155 * integration)
156 * @param vecAbsoluteTolerance allowed absolute error
157 * @param vecRelativeTolerance allowed relative error
158 */
159 protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal,
160 final double minStep, final double maxStep,
161 final double[] vecAbsoluteTolerance,
162 final double[] vecRelativeTolerance) {
163
164 super(field, name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
165 this.usingFieldCoefficients = false;
166
167 this.fsal = fsal;
168 this.c = getC();
169 this.a = getA();
170 this.b = getB();
171
172 exp = -1.0 / getOrder();
173
174 // set the default values of the algorithm control parameters
175 setSafety(field.getZero().add(0.9));
176 setMinReduction(field.getZero().add(0.2));
177 setMaxGrowth(field.getZero().add(10.0));
178
179 }
180
181 /** Create an interpolator.
182 * @param forward integration direction indicator
183 * @param yDotK slopes at the intermediate points
184 * @param globalPreviousState start of the global step
185 * @param globalCurrentState end of the global step
186 * @param mapper equations mapper for the all equations
187 * @return external weights for the high order method from Butcher array
188 */
189 protected abstract RungeKuttaFieldStateInterpolator<T> createInterpolator(boolean forward, T[][] yDotK,
190 FieldODEStateAndDerivative<T> globalPreviousState,
191 FieldODEStateAndDerivative<T> globalCurrentState,
192 FieldEquationsMapper<T> mapper);
193
194 /** Get the order of the method.
195 * @return order of the method
196 */
197 public abstract int getOrder();
198
199 /** Get the safety factor for stepsize control.
200 * @return safety factor
201 */
202 public T getSafety() {
203 return safety;
204 }
205
206 /** Set the safety factor for stepsize control.
207 * @param safety safety factor
208 */
209 public void setSafety(final T safety) {
210 this.safety = safety;
211 }
212
213 /**
214 * Setter for the flag between real or Field coefficients in the Butcher array.
215 *
216 * @param usingFieldCoefficients new value for flag
217 */
218 public void setUsingFieldCoefficients(boolean usingFieldCoefficients) {
219 this.usingFieldCoefficients = usingFieldCoefficients;
220 }
221
222 /** {@inheritDoc} */
223 @Override
224 public boolean isUsingFieldCoefficients() {
225 return usingFieldCoefficients;
226 }
227
228 /** {@inheritDoc} */
229 @Override
230 public int getNumberOfStages() {
231 return b.length;
232 }
233
234 /** {@inheritDoc} */
235 @Override
236 protected FieldODEStateAndDerivative<T> initIntegration(FieldExpandableODE<T> eqn, FieldODEState<T> s0, T t) {
237 if (!isUsingFieldCoefficients()) {
238 realA = getRealA();
239 realB = getRealB();
240 realC = getRealC();
241 }
242 return super.initIntegration(eqn, s0, t);
243 }
244
245 /** {@inheritDoc} */
246 @Override
247 public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
248 final FieldODEState<T> initialState, final T finalTime)
249 throws MathIllegalArgumentException, MathIllegalStateException {
250
251 sanityChecks(initialState, finalTime);
252 setStepStart(initIntegration(equations, initialState, finalTime));
253 final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0;
254
255 // create some internal working arrays
256 final int stages = getNumberOfStages();
257 final T[][] yDotK = MathArrays.buildArray(getField(), stages, -1);
258 T[] yTmp = MathArrays.buildArray(getField(), equations.getMapper().getTotalDimension());
259
260 // set up integration control objects
261 T hNew = getField().getZero();
262 boolean firstTime = true;
263
264 // main integration loop
265 setIsLastStep(false);
266 do {
267
268 // iterate over step size, ensuring local normalized error is smaller than 1
269 double error = 10.0;
270 while (error >= 1.0) {
271
272 // first stage
273 final T[] y = getStepStart().getCompleteState();
274 yDotK[0] = getStepStart().getCompleteDerivative();
275
276 if (firstTime) {
277 final StepsizeHelper helper = getStepSizeHelper();
278 final T[] scale = MathArrays.buildArray(getField(), helper.getMainSetDimension());
279 for (int i = 0; i < scale.length; ++i) {
280 scale[i] = helper.getTolerance(i, y[i].abs());
281 }
282 hNew = getField().getZero().add(initializeStep(forward, getOrder(), scale, getStepStart(), equations.getMapper()));
283 firstTime = false;
284 }
285
286 setStepSize(hNew);
287 if (forward) {
288 if (getStepStart().getTime().add(getStepSize()).subtract(finalTime).getReal() >= 0) {
289 setStepSize(finalTime.subtract(getStepStart().getTime()));
290 }
291 } else {
292 if (getStepStart().getTime().add(getStepSize()).subtract(finalTime).getReal() <= 0) {
293 setStepSize(finalTime.subtract(getStepStart().getTime()));
294 }
295 }
296
297 // next stages
298 if (isUsingFieldCoefficients()) {
299 FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(),
300 getStepStart().getTime(), y, getStepSize(), a, c, yDotK);
301 yTmp = FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y, yDotK, getStepSize(), b);
302 } else {
303 FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(),
304 getStepStart().getTime(), y, getStepSize(), realA, realC, yDotK);
305 yTmp = FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y, yDotK, getStepSize(), realB);
306 }
307
308 incrementEvaluations(stages - 1);
309
310 // estimate the error at the end of the step
311 error = estimateError(yDotK, y, yTmp, getStepSize());
312 if (error >= 1.0) {
313 // reject the step and attempt to reduce error by stepsize control
314 final T factor = MathUtils.min(maxGrowth,
315 MathUtils.max(minReduction, safety.multiply(FastMath.pow(error, exp))));
316 hNew = getStepSizeHelper().filterStep(getStepSize().multiply(factor), forward, false);
317 }
318
319 }
320 final T stepEnd = getStepStart().getTime().add(getStepSize());
321 final T[] yDotTmp = (fsal >= 0) ? yDotK[fsal] : computeDerivatives(stepEnd, yTmp);
322 final FieldODEStateAndDerivative<T> stateTmp = equations.getMapper().mapStateAndDerivative(stepEnd, yTmp, yDotTmp);
323
324 // local error is small enough: accept the step, trigger events and step handlers
325 setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp, equations.getMapper()),
326 finalTime));
327
328 if (!isLastStep()) {
329
330 // stepsize control for next step
331 final T factor = MathUtils.min(maxGrowth,
332 MathUtils.max(minReduction, safety.multiply(FastMath.pow(error, exp))));
333 final T scaledH = getStepSize().multiply(factor);
334 final T nextT = getStepStart().getTime().add(scaledH);
335 final boolean nextIsLast = forward ?
336 nextT.subtract(finalTime).getReal() >= 0 :
337 nextT.subtract(finalTime).getReal() <= 0;
338 hNew = getStepSizeHelper().filterStep(scaledH, forward, nextIsLast);
339
340 final T filteredNextT = getStepStart().getTime().add(hNew);
341 final boolean filteredNextIsLast = forward ?
342 filteredNextT.subtract(finalTime).getReal() >= 0 :
343 filteredNextT.subtract(finalTime).getReal() <= 0;
344 if (filteredNextIsLast) {
345 hNew = finalTime.subtract(getStepStart().getTime());
346 }
347
348 }
349
350 } while (!isLastStep());
351
352 final FieldODEStateAndDerivative<T> finalState = getStepStart();
353 resetInternalState();
354 return finalState;
355
356 }
357
358 /** Get the minimal reduction factor for stepsize control.
359 * @return minimal reduction factor
360 */
361 public T getMinReduction() {
362 return minReduction;
363 }
364
365 /** Set the minimal reduction factor for stepsize control.
366 * @param minReduction minimal reduction factor
367 */
368 public void setMinReduction(final T minReduction) {
369 this.minReduction = minReduction;
370 }
371
372 /** Get the maximal growth factor for stepsize control.
373 * @return maximal growth factor
374 */
375 public T getMaxGrowth() {
376 return maxGrowth;
377 }
378
379 /** Set the maximal growth factor for stepsize control.
380 * @param maxGrowth maximal growth factor
381 */
382 public void setMaxGrowth(final T maxGrowth) {
383 this.maxGrowth = maxGrowth;
384 }
385
386 /** Compute the error ratio.
387 * @param yDotK derivatives computed during the first stages
388 * @param y0 estimate of the step at the start of the step
389 * @param y1 estimate of the step at the end of the step
390 * @param h current step
391 * @return error ratio, greater than 1 if step should be rejected
392 */
393 protected abstract double estimateError(T[][] yDotK, T[] y0, T[] y1, T h);
394
395 }