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
21 import org.hipparchus.CalculusFieldElement;
22 import org.hipparchus.Field;
23 import org.hipparchus.ode.FieldExpandableODE;
24 import org.hipparchus.ode.FieldODEIntegrator;
25 import org.hipparchus.ode.FieldOrdinaryDifferentialEquation;
26 import org.hipparchus.util.MathArrays;
27
28 /**
29 * This interface implements the part of Runge-Kutta
30 * Field integrators for Ordinary Differential Equations
31 * common to fixed- and adaptive steps.
32 *
33 * <p>These methods are explicit Runge-Kutta methods, their Butcher
34 * arrays are as follows :</p>
35 * <pre>
36 * 0 |
37 * c2 | a21
38 * c3 | a31 a32
39 * ... | ...
40 * cs | as1 as2 ... ass-1
41 * |--------------------------
42 * | b1 b2 ... bs-1 bs
43 * </pre>
44 *
45 * @see FieldButcherArrayProvider
46 * @see FixedStepRungeKuttaFieldIntegrator
47 * @see EmbeddedRungeKuttaFieldIntegrator
48 * @param <T> the type of the field elements
49 * @since 3.1
50 */
51
52 public interface FieldExplicitRungeKuttaIntegrator<T extends CalculusFieldElement<T>>
53 extends FieldButcherArrayProvider<T>, FieldODEIntegrator<T> {
54
55 /** Get the time steps from Butcher array (without the first zero). Real version (non-Field).
56 * @return time steps from Butcher array (without the first zero).
57 */
58 default double[] getRealC() {
59 final T[] c = getC();
60 final double[] cReal = new double[c.length];
61 for (int i = 0; i < c.length; i++) {
62 cReal[i] = c[i].getReal();
63 }
64 return cReal;
65 }
66
67 /** Get the internal weights from Butcher array (without the first empty row). Real version (non-Field).
68 * @return internal weights from Butcher array (without the first empty row)
69 */
70 default double[][] getRealA() {
71 final T[][] a = getA();
72 final double[][] aReal = new double[a.length][];
73 for (int i = 0; i < a.length; i++) {
74 aReal[i] = new double[a[i].length];
75 for (int j = 0; j < aReal[i].length; j++) {
76 aReal[i][j] = a[i][j].getReal();
77 }
78 }
79 return aReal;
80 }
81
82 /** Get the external weights for the high order method from Butcher array. Real version (non-Field).
83 * @return external weights for the high order method from Butcher array
84 */
85 default double[] getRealB() {
86 final T[] b = getB();
87 final double[] bReal = new double[b.length];
88 for (int i = 0; i < b.length; i++) {
89 bReal[i] = b[i].getReal();
90 }
91 return bReal;
92 }
93
94 /**
95 * Getter for the flag between real or Field coefficients in the Butcher array.
96 *
97 * @return flag
98 */
99 boolean isUsingFieldCoefficients();
100
101 /**
102 * Getter for the number of stages corresponding to the Butcher array.
103 *
104 * @return number of stages
105 */
106 default int getNumberOfStages() {
107 return getB().length;
108 }
109
110 /** Fast computation of a single step of ODE integration.
111 * <p>This method is intended for the limited use case of
112 * very fast computation of only one step without using any of the
113 * rich features of general integrators that may take some time
114 * to set up (i.e. no step handlers, no events handlers, no additional
115 * states, no interpolators, no error control, no evaluations count,
116 * no sanity checks ...). It handles the strict minimum of computation,
117 * so it can be embedded in outer loops.</p>
118 * <p>
119 * This method is <em>not</em> used at all by the {@link #integrate(FieldExpandableODE,
120 * org.hipparchus.ode.FieldODEState, CalculusFieldElement)} method. It also completely ignores the step set at
121 * construction time, and uses only a single step to go from {@code t0} to {@code t}.
122 * </p>
123 * <p>
124 * As this method does not use any of the state-dependent features of the integrator,
125 * it should be reasonably thread-safe <em>if and only if</em> the provided differential
126 * equations are themselves thread-safe.
127 * </p>
128 * @param equations differential equations to integrate
129 * @param t0 initial time
130 * @param y0 initial value of the state vector at t0
131 * @param t target time for the integration
132 * (can be set to a value smaller than {@code t0} for backward integration)
133 * @return state vector at {@code t}
134 */
135 default T[] singleStep(final FieldOrdinaryDifferentialEquation<T> equations, final T t0, final T[] y0, final T t) {
136
137 // create some internal working arrays
138 final int stages = getNumberOfStages();
139 final T[][] yDotK = MathArrays.buildArray(t0.getField(), stages, -1);
140
141 // first stage
142 final T h = t.subtract(t0);
143 final FieldExpandableODE<T> fieldExpandableODE = new FieldExpandableODE<>(equations);
144 yDotK[0] = fieldExpandableODE.computeDerivatives(t0, y0);
145
146 if (isUsingFieldCoefficients()) {
147 applyInternalButcherWeights(fieldExpandableODE, t0, y0, h, getA(), getC(), yDotK);
148 return applyExternalButcherWeights(y0, yDotK, h, getB());
149 } else {
150 applyInternalButcherWeights(fieldExpandableODE, t0, y0, h, getRealA(), getRealC(), yDotK);
151 return applyExternalButcherWeights(y0, yDotK, h, getRealB());
152 }
153 }
154
155 /**
156 * Create a fraction from integers.
157 *
158 * @param <T> the type of the field elements
159 * @param field field to which elements belong
160 * @param p numerator
161 * @param q denominator
162 * @return p/q computed in the instance field
163 */
164 static <T extends CalculusFieldElement<T>> T fraction(final Field<T> field, final int p, final int q) {
165 final T zero = field.getZero();
166 return zero.newInstance(p).divide(zero.newInstance(q));
167 }
168
169 /**
170 * Create a fraction from doubles.
171 * @param <T> the type of the field elements
172 * @param field field to which elements belong
173 * @param p numerator
174 * @param q denominator
175 * @return p/q computed in the instance field
176 */
177 static <T extends CalculusFieldElement<T>> T fraction(final Field<T> field, final double p, final double q) {
178 final T zero = field.getZero();
179 return zero.newInstance(p).divide(zero.newInstance(q));
180 }
181
182 /**
183 * Apply internal weights of Butcher array, with corresponding times.
184 * @param <T> the type of the field elements
185 * @param equations differential equations to integrate
186 * @param t0 initial time
187 * @param y0 initial value of the state vector at t0
188 * @param h step size
189 * @param a internal weights of Butcher array
190 * @param c times of Butcher array
191 * @param yDotK array where to store result
192 */
193 static <T extends CalculusFieldElement<T>> void applyInternalButcherWeights(final FieldExpandableODE<T> equations,
194 final T t0, final T[] y0, final T h,
195 final T[][] a, final T[] c,
196 final T[][] yDotK) {
197 // create some internal working arrays
198 final int stages = c.length + 1;
199 final T[] yTmp = y0.clone();
200
201 for (int k = 1; k < stages; ++k) {
202
203 for (int j = 0; j < y0.length; ++j) {
204 T sum = yDotK[0][j].multiply(a[k - 1][0]);
205 for (int l = 1; l < k; ++l) {
206 sum = sum.add(yDotK[l][j].multiply(a[k - 1][l]));
207 }
208 yTmp[j] = y0[j].add(h.multiply(sum));
209 }
210
211 yDotK[k] = equations.computeDerivatives(t0.add(h.multiply(c[k - 1])), yTmp);
212 }
213 }
214
215 /** Apply internal weights of Butcher array, with corresponding times. Version with real Butcher array (non-Field).
216 * @param <T> the type of the field elements
217 * @param equations differential equations to integrate
218 * @param t0 initial time
219 * @param y0 initial value of the state vector at t0
220 * @param h step size
221 * @param a internal weights of Butcher array
222 * @param c times of Butcher array
223 * @param yDotK array where to store result
224 */
225 static <T extends CalculusFieldElement<T>> void applyInternalButcherWeights(final FieldExpandableODE<T> equations,
226 final T t0, final T[] y0, final T h,
227 final double[][] a, final double[] c,
228 final T[][] yDotK) {
229 // create some internal working arrays
230 final int stages = c.length + 1;
231 final T[] yTmp = y0.clone();
232
233 for (int k = 1; k < stages; ++k) {
234
235 for (int j = 0; j < y0.length; ++j) {
236 T sum = yDotK[0][j].multiply(a[k - 1][0]);
237 for (int l = 1; l < k; ++l) {
238 sum = sum.add(yDotK[l][j].multiply(a[k - 1][l]));
239 }
240 yTmp[j] = y0[j].add(h.multiply(sum));
241 }
242
243 yDotK[k] = equations.computeDerivatives(t0.add(h.multiply(c[k - 1])), yTmp);
244 }
245 }
246
247 /** Apply external weights of Butcher array, assuming internal ones have been applied.
248 * @param <T> the type of the field elements
249 * @param yDotK output of stages
250 * @param y0 initial value of the state vector at t0
251 * @param h step size
252 * @param b external weights of Butcher array
253 * @return state vector
254 */
255 static <T extends CalculusFieldElement<T>> T[] applyExternalButcherWeights(final T[] y0, final T[][] yDotK,
256 final T h, final T[] b) {
257 final T[] y = y0.clone();
258 final int stages = b.length;
259 for (int j = 0; j < y0.length; ++j) {
260 T sum = yDotK[0][j].multiply(b[0]);
261 for (int l = 1; l < stages; ++l) {
262 sum = sum.add(yDotK[l][j].multiply(b[l]));
263 }
264 y[j] = y[j].add(h.multiply(sum));
265 }
266 return y;
267 }
268
269 /** Apply external weights of Butcher array, assuming internal ones have been applied. Version with real Butcher
270 * array (non-Field version).
271 * @param <T> the type of the field elements
272 * @param yDotK output of stages
273 * @param y0 initial value of the state vector at t0
274 * @param h step size
275 * @param b external weights of Butcher array
276 * @return state vector
277 */
278 static <T extends CalculusFieldElement<T>> T[] applyExternalButcherWeights(final T[] y0, final T[][] yDotK,
279 final T h, final double[] b) {
280 final T[] y = y0.clone();
281 final int stages = b.length;
282 for (int j = 0; j < y0.length; ++j) {
283 T sum = yDotK[0][j].multiply(b[0]);
284 for (int l = 1; l < stages; ++l) {
285 sum = sum.add(yDotK[l][j].multiply(b[l]));
286 }
287 y[j] = y[j].add(h.multiply(sum));
288 }
289 return y;
290 }
291
292 }