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