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