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
21 import org.apache.commons.math4.legacy.core.Field;
22 import org.apache.commons.math4.legacy.core.RealFieldElement;
23 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
24 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
25 import org.apache.commons.math4.legacy.exception.NoBracketingException;
26 import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
27 import org.apache.commons.math4.legacy.ode.AbstractFieldIntegrator;
28 import org.apache.commons.math4.legacy.ode.FieldEquationsMapper;
29 import org.apache.commons.math4.legacy.ode.FieldExpandableODE;
30 import org.apache.commons.math4.legacy.ode.FirstOrderFieldDifferentialEquations;
31 import org.apache.commons.math4.legacy.ode.FieldODEState;
32 import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative;
33 import org.apache.commons.math4.legacy.core.MathArrays;
34
35 /**
36 * This class implements the common part of all fixed step Runge-Kutta
37 * integrators for Ordinary Differential Equations.
38 *
39 * <p>These methods are explicit Runge-Kutta methods, their Butcher
40 * arrays are as follows :
41 * <pre>
42 * 0 |
43 * c2 | a21
44 * c3 | a31 a32
45 * ... | ...
46 * cs | as1 as2 ... ass-1
47 * |--------------------------
48 * | b1 b2 ... bs-1 bs
49 * </pre>
50 *
51 * @see EulerFieldIntegrator
52 * @see ClassicalRungeKuttaFieldIntegrator
53 * @see GillFieldIntegrator
54 * @see MidpointFieldIntegrator
55 * @param <T> the type of the field elements
56 * @since 3.6
57 */
58
59 public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
60 extends AbstractFieldIntegrator<T>
61 implements FieldButcherArrayProvider<T> {
62
63 /** Time steps from Butcher array (without the first zero). */
64 private final T[] c;
65
66 /** Internal weights from Butcher array (without the first empty row). */
67 private final T[][] a;
68
69 /** External weights for the high order method from Butcher array. */
70 private final T[] b;
71
72 /** Integration step. */
73 private final T step;
74
75 /** Simple constructor.
76 * Build a Runge-Kutta integrator with the given
77 * step. The default step handler does nothing.
78 * @param field field to which the time and state vector elements belong
79 * @param name name of the method
80 * @param step integration step
81 */
82 protected RungeKuttaFieldIntegrator(final Field<T> field, final String name, final T step) {
83 super(field, name);
84 this.c = getC();
85 this.a = getA();
86 this.b = getB();
87 this.step = step.abs();
88 }
89
90 /** Create a fraction.
91 * @param p numerator
92 * @param q denominator
93 * @return p/q computed in the instance field
94 */
95 protected T fraction(final int p, final int q) {
96 return getField().getZero().add(p).divide(q);
97 }
98
99 /** Create an interpolator.
100 * @param forward integration direction indicator
101 * @param yDotK slopes at the intermediate points
102 * @param globalPreviousState start of the global step
103 * @param globalCurrentState end of the global step
104 * @param mapper equations mapper for the all equations
105 * @return external weights for the high order method from Butcher array
106 */
107 protected abstract RungeKuttaFieldStepInterpolator<T> createInterpolator(boolean forward, T[][] yDotK,
108 FieldODEStateAndDerivative<T> globalPreviousState,
109 FieldODEStateAndDerivative<T> globalCurrentState,
110 FieldEquationsMapper<T> mapper);
111
112 /** {@inheritDoc} */
113 @Override
114 public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
115 final FieldODEState<T> initialState, final T finalTime)
116 throws NumberIsTooSmallException, DimensionMismatchException,
117 MaxCountExceededException, NoBracketingException {
118
119 sanityChecks(initialState, finalTime);
120 final T t0 = initialState.getTime();
121 final T[] y0 = equations.getMapper().mapState(initialState);
122 setStepStart(initIntegration(equations, t0, y0, finalTime));
123 final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0;
124
125 // create some internal working arrays
126 final int stages = c.length + 1;
127 T[] y = y0;
128 final T[][] yDotK = MathArrays.buildArray(getField(), stages, -1);
129 final T[] yTmp = MathArrays.buildArray(getField(), y0.length);
130
131 // set up integration control objects
132 if (forward) {
133 if (getStepStart().getTime().add(step).subtract(finalTime).getReal() >= 0) {
134 setStepSize(finalTime.subtract(getStepStart().getTime()));
135 } else {
136 setStepSize(step);
137 }
138 } else {
139 if (getStepStart().getTime().subtract(step).subtract(finalTime).getReal() <= 0) {
140 setStepSize(finalTime.subtract(getStepStart().getTime()));
141 } else {
142 setStepSize(step.negate());
143 }
144 }
145
146 // main integration loop
147 setIsLastStep(false);
148 do {
149
150 // first stage
151 y = equations.getMapper().mapState(getStepStart());
152 yDotK[0] = equations.getMapper().mapDerivative(getStepStart());
153
154 // next stages
155 for (int k = 1; k < stages; ++k) {
156
157 for (int j = 0; j < y0.length; ++j) {
158 T sum = yDotK[0][j].multiply(a[k-1][0]);
159 for (int l = 1; l < k; ++l) {
160 sum = sum.add(yDotK[l][j].multiply(a[k-1][l]));
161 }
162 yTmp[j] = y[j].add(getStepSize().multiply(sum));
163 }
164
165 yDotK[k] = computeDerivatives(getStepStart().getTime().add(getStepSize().multiply(c[k-1])), yTmp);
166 }
167
168 // estimate the state at the end of the step
169 for (int j = 0; j < y0.length; ++j) {
170 T sum = yDotK[0][j].multiply(b[0]);
171 for (int l = 1; l < stages; ++l) {
172 sum = sum.add(yDotK[l][j].multiply(b[l]));
173 }
174 yTmp[j] = y[j].add(getStepSize().multiply(sum));
175 }
176 final T stepEnd = getStepStart().getTime().add(getStepSize());
177 final T[] yDotTmp = computeDerivatives(stepEnd, yTmp);
178 final FieldODEStateAndDerivative<T> stateTmp = new FieldODEStateAndDerivative<>(stepEnd, yTmp, yDotTmp);
179
180 // discrete events handling
181 System.arraycopy(yTmp, 0, y, 0, y0.length);
182 setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp, equations.getMapper()),
183 finalTime));
184
185 if (!isLastStep()) {
186
187 // stepsize control for next step
188 final T nextT = getStepStart().getTime().add(getStepSize());
189 final boolean nextIsLast = forward ?
190 (nextT.subtract(finalTime).getReal() >= 0) :
191 (nextT.subtract(finalTime).getReal() <= 0);
192 if (nextIsLast) {
193 setStepSize(finalTime.subtract(getStepStart().getTime()));
194 }
195 }
196 } while (!isLastStep());
197
198 final FieldODEStateAndDerivative<T> finalState = getStepStart();
199 setStepStart(null);
200 setStepSize(null);
201 return finalState;
202 }
203
204 /** Fast computation of a single step of ODE integration.
205 * <p>This method is intended for the limited use case of
206 * very fast computation of only one step without using any of the
207 * rich features of general integrators that may take some time
208 * to set up (i.e. no step handlers, no events handlers, no additional
209 * states, no interpolators, no error control, no evaluations count,
210 * no sanity checks ...). It handles the strict minimum of computation,
211 * so it can be embedded in outer loops.</p>
212 * <p>
213 * This method is <em>not</em> used at all by the {@link #integrate(FieldExpandableODE,
214 * FieldODEState, RealFieldElement)} method. It also completely ignores the step set at
215 * construction time, and uses only a single step to go from {@code t0} to {@code t}.
216 * </p>
217 * <p>
218 * As this method does not use any of the state-dependent features of the integrator,
219 * it should be reasonably thread-safe <em>if and only if</em> the provided differential
220 * equations are themselves thread-safe.
221 * </p>
222 * @param equations differential equations to integrate
223 * @param t0 initial time
224 * @param y0 initial value of the state vector at t0
225 * @param t target time for the integration
226 * (can be set to a value smaller than {@code t0} for backward integration)
227 * @return state vector at {@code t}
228 */
229 public T[] singleStep(final FirstOrderFieldDifferentialEquations<T> equations,
230 final T t0, final T[] y0, final T t) {
231
232 // create some internal working arrays
233 final T[] y = y0.clone();
234 final int stages = c.length + 1;
235 final T[][] yDotK = MathArrays.buildArray(getField(), stages, -1);
236 final T[] yTmp = y0.clone();
237
238 // first stage
239 final T h = t.subtract(t0);
240 yDotK[0] = equations.computeDerivatives(t0, y);
241
242 // next stages
243 for (int k = 1; k < stages; ++k) {
244
245 for (int j = 0; j < y0.length; ++j) {
246 T sum = yDotK[0][j].multiply(a[k-1][0]);
247 for (int l = 1; l < k; ++l) {
248 sum = sum.add(yDotK[l][j].multiply(a[k-1][l]));
249 }
250 yTmp[j] = y[j].add(h.multiply(sum));
251 }
252
253 yDotK[k] = equations.computeDerivatives(t0.add(h.multiply(c[k-1])), yTmp);
254 }
255
256 // estimate the state at the end of the step
257 for (int j = 0; j < y0.length; ++j) {
258 T sum = yDotK[0][j].multiply(b[0]);
259 for (int l = 1; l < stages; ++l) {
260 sum = sum.add(yDotK[l][j].multiply(b[l]));
261 }
262 y[j] = y[j].add(h.multiply(sum));
263 }
264
265 return y;
266 }
267 }