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 RungeKuttaFieldIntegrator
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 FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(fieldExpandableODE, t0, y0, h, getA(), getC(),
148 yDotK);
149 return FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y0, yDotK, h, getB());
150 } else {
151 FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(fieldExpandableODE, t0, y0, h,
152 getRealA(), getRealC(), yDotK);
153 return FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y0, yDotK, h, getRealB());
154 }
155 }
156
157 /**
158 * Create a fraction from integers.
159 *
160 * @param <T> the type of the field elements
161 * @param field field to which elements belong
162 * @param p numerator
163 * @param q denominator
164 * @return p/q computed in the instance field
165 */
166 static <T extends CalculusFieldElement<T>> T fraction(final Field<T> field, final int p, final int q) {
167 final T zero = field.getZero();
168 return zero.newInstance(p).divide(zero.newInstance(q));
169 }
170
171 /**
172 * Create a fraction from doubles.
173 * @param <T> the type of the field elements
174 * @param field field to which elements belong
175 * @param p numerator
176 * @param q denominator
177 * @return p/q computed in the instance field
178 */
179 static <T extends CalculusFieldElement<T>> T fraction(final Field<T> field, final double p, final double q) {
180 final T zero = field.getZero();
181 return zero.newInstance(p).divide(zero.newInstance(q));
182 }
183
184 /**
185 * Apply internal weights of Butcher array, with corresponding times.
186 * @param <T> the type of the field elements
187 * @param equations differential equations to integrate
188 * @param t0 initial time
189 * @param y0 initial value of the state vector at t0
190 * @param h step size
191 * @param a internal weights of Butcher array
192 * @param c times of Butcher array
193 * @param yDotK array where to store result
194 */
195 static <T extends CalculusFieldElement<T>> void applyInternalButcherWeights(final FieldExpandableODE<T> equations,
196 final T t0, final T[] y0, final T h,
197 final T[][] a, final T[] c,
198 final T[][] yDotK) {
199 // create some internal working arrays
200 final int stages = c.length + 1;
201 final T[] yTmp = y0.clone();
202
203 for (int k = 1; k < stages; ++k) {
204
205 for (int j = 0; j < y0.length; ++j) {
206 T sum = yDotK[0][j].multiply(a[k - 1][0]);
207 for (int l = 1; l < k; ++l) {
208 sum = sum.add(yDotK[l][j].multiply(a[k - 1][l]));
209 }
210 yTmp[j] = y0[j].add(h.multiply(sum));
211 }
212
213 yDotK[k] = equations.computeDerivatives(t0.add(h.multiply(c[k - 1])), yTmp);
214 }
215 }
216
217 /** Apply internal weights of Butcher array, with corresponding times. Version with real Butcher array (non-Field).
218 * @param <T> the type of the field elements
219 * @param equations differential equations to integrate
220 * @param t0 initial time
221 * @param y0 initial value of the state vector at t0
222 * @param h step size
223 * @param a internal weights of Butcher array
224 * @param c times of Butcher array
225 * @param yDotK array where to store result
226 */
227 static <T extends CalculusFieldElement<T>> void applyInternalButcherWeights(final FieldExpandableODE<T> equations,
228 final T t0, final T[] y0, final T h,
229 final double[][] a, final double[] c,
230 final T[][] yDotK) {
231 // create some internal working arrays
232 final int stages = c.length + 1;
233 final T[] yTmp = y0.clone();
234
235 for (int k = 1; k < stages; ++k) {
236
237 for (int j = 0; j < y0.length; ++j) {
238 T sum = yDotK[0][j].multiply(a[k - 1][0]);
239 for (int l = 1; l < k; ++l) {
240 sum = sum.add(yDotK[l][j].multiply(a[k - 1][l]));
241 }
242 yTmp[j] = y0[j].add(h.multiply(sum));
243 }
244
245 yDotK[k] = equations.computeDerivatives(t0.add(h.multiply(c[k - 1])), yTmp);
246 }
247 }
248
249 /** Apply external weights of Butcher array, assuming internal ones have been applied.
250 * @param <T> the type of the field elements
251 * @param yDotK output of stages
252 * @param y0 initial value of the state vector at t0
253 * @param h step size
254 * @param b external weights of Butcher array
255 * @return state vector
256 */
257 static <T extends CalculusFieldElement<T>> T[] applyExternalButcherWeights(final T[] y0, final T[][] yDotK,
258 final T h, final T[] b) {
259 final T[] y = y0.clone();
260 final int stages = b.length;
261 for (int j = 0; j < y0.length; ++j) {
262 T sum = yDotK[0][j].multiply(b[0]);
263 for (int l = 1; l < stages; ++l) {
264 sum = sum.add(yDotK[l][j].multiply(b[l]));
265 }
266 y[j] = y[j].add(h.multiply(sum));
267 }
268 return y;
269 }
270
271 /** Apply external weights of Butcher array, assuming internal ones have been applied. Version with real Butcher
272 * array (non-Field version).
273 * @param <T> the type of the field elements
274 * @param yDotK output of stages
275 * @param y0 initial value of the state vector at t0
276 * @param h step size
277 * @param b external weights of Butcher array
278 * @return state vector
279 */
280 static <T extends CalculusFieldElement<T>> T[] applyExternalButcherWeights(final T[] y0, final T[][] yDotK,
281 final T h, final double[] b) {
282 final T[] y = y0.clone();
283 final int stages = b.length;
284 for (int j = 0; j < y0.length; ++j) {
285 T sum = yDotK[0][j].multiply(b[0]);
286 for (int l = 1; l < stages; ++l) {
287 sum = sum.add(yDotK[l][j].multiply(b[l]));
288 }
289 y[j] = y[j].add(h.multiply(sum));
290 }
291 return y;
292 }
293
294 }