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