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