FieldExplicitRungeKuttaIntegrator.java
/*
* Licensed to the Hipparchus project under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The Hipparchus project licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* https://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.hipparchus.ode.nonstiff;
import org.hipparchus.CalculusFieldElement;
import org.hipparchus.Field;
import org.hipparchus.ode.FieldExpandableODE;
import org.hipparchus.ode.FieldODEIntegrator;
import org.hipparchus.ode.FieldOrdinaryDifferentialEquation;
import org.hipparchus.util.MathArrays;
/**
* This interface implements the part of Runge-Kutta
* Field integrators for Ordinary Differential Equations
* common to fixed- and adaptive steps.
*
* <p>These methods are explicit Runge-Kutta methods, their Butcher
* arrays are as follows :</p>
* <pre>
* 0 |
* c2 | a21
* c3 | a31 a32
* ... | ...
* cs | as1 as2 ... ass-1
* |--------------------------
* | b1 b2 ... bs-1 bs
* </pre>
*
* @see FieldButcherArrayProvider
* @see RungeKuttaFieldIntegrator
* @see EmbeddedRungeKuttaFieldIntegrator
* @param <T> the type of the field elements
* @since 3.1
*/
public interface FieldExplicitRungeKuttaIntegrator<T extends CalculusFieldElement<T>>
extends FieldButcherArrayProvider<T>, FieldODEIntegrator<T> {
/** Get the time steps from Butcher array (without the first zero). Real version (non-Field).
* @return time steps from Butcher array (without the first zero).
*/
default double[] getRealC() {
final T[] c = getC();
final double[] cReal = new double[c.length];
for (int i = 0; i < c.length; i++) {
cReal[i] = c[i].getReal();
}
return cReal;
}
/** Get the internal weights from Butcher array (without the first empty row). Real version (non-Field).
* @return internal weights from Butcher array (without the first empty row)
*/
default double[][] getRealA() {
final T[][] a = getA();
final double[][] aReal = new double[a.length][];
for (int i = 0; i < a.length; i++) {
aReal[i] = new double[a[i].length];
for (int j = 0; j < aReal[i].length; j++) {
aReal[i][j] = a[i][j].getReal();
}
}
return aReal;
}
/** Get the external weights for the high order method from Butcher array. Real version (non-Field).
* @return external weights for the high order method from Butcher array
*/
default double[] getRealB() {
final T[] b = getB();
final double[] bReal = new double[b.length];
for (int i = 0; i < b.length; i++) {
bReal[i] = b[i].getReal();
}
return bReal;
}
/**
* Getter for the flag between real or Field coefficients in the Butcher array.
*
* @return flag
*/
boolean isUsingFieldCoefficients();
/**
* Getter for the number of stages corresponding to the Butcher array.
*
* @return number of stages
*/
default int getNumberOfStages() {
return getB().length;
}
/** Fast computation of a single step of ODE integration.
* <p>This method is intended for the limited use case of
* very fast computation of only one step without using any of the
* rich features of general integrators that may take some time
* to set up (i.e. no step handlers, no events handlers, no additional
* states, no interpolators, no error control, no evaluations count,
* no sanity checks ...). It handles the strict minimum of computation,
* so it can be embedded in outer loops.</p>
* <p>
* This method is <em>not</em> used at all by the {@link #integrate(FieldExpandableODE,
* org.hipparchus.ode.FieldODEState, CalculusFieldElement)} method. It also completely ignores the step set at
* construction time, and uses only a single step to go from {@code t0} to {@code t}.
* </p>
* <p>
* As this method does not use any of the state-dependent features of the integrator,
* it should be reasonably thread-safe <em>if and only if</em> the provided differential
* equations are themselves thread-safe.
* </p>
* @param equations differential equations to integrate
* @param t0 initial time
* @param y0 initial value of the state vector at t0
* @param t target time for the integration
* (can be set to a value smaller than {@code t0} for backward integration)
* @return state vector at {@code t}
*/
default T[] singleStep(final FieldOrdinaryDifferentialEquation<T> equations, final T t0, final T[] y0, final T t) {
// create some internal working arrays
final int stages = getNumberOfStages();
final T[][] yDotK = MathArrays.buildArray(t0.getField(), stages, -1);
// first stage
final T h = t.subtract(t0);
final FieldExpandableODE<T> fieldExpandableODE = new FieldExpandableODE<>(equations);
yDotK[0] = fieldExpandableODE.computeDerivatives(t0, y0);
if (isUsingFieldCoefficients()) {
FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(fieldExpandableODE, t0, y0, h, getA(), getC(),
yDotK);
return FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y0, yDotK, h, getB());
} else {
FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(fieldExpandableODE, t0, y0, h,
getRealA(), getRealC(), yDotK);
return FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y0, yDotK, h, getRealB());
}
}
/**
* Create a fraction from integers.
*
* @param <T> the type of the field elements
* @param field field to which elements belong
* @param p numerator
* @param q denominator
* @return p/q computed in the instance field
*/
static <T extends CalculusFieldElement<T>> T fraction(final Field<T> field, final int p, final int q) {
final T zero = field.getZero();
return zero.newInstance(p).divide(zero.newInstance(q));
}
/**
* Create a fraction from doubles.
* @param <T> the type of the field elements
* @param field field to which elements belong
* @param p numerator
* @param q denominator
* @return p/q computed in the instance field
*/
static <T extends CalculusFieldElement<T>> T fraction(final Field<T> field, final double p, final double q) {
final T zero = field.getZero();
return zero.newInstance(p).divide(zero.newInstance(q));
}
/**
* Apply internal weights of Butcher array, with corresponding times.
* @param <T> the type of the field elements
* @param equations differential equations to integrate
* @param t0 initial time
* @param y0 initial value of the state vector at t0
* @param h step size
* @param a internal weights of Butcher array
* @param c times of Butcher array
* @param yDotK array where to store result
*/
static <T extends CalculusFieldElement<T>> void applyInternalButcherWeights(final FieldExpandableODE<T> equations,
final T t0, final T[] y0, final T h,
final T[][] a, final T[] c,
final T[][] yDotK) {
// create some internal working arrays
final int stages = c.length + 1;
final T[] yTmp = y0.clone();
for (int k = 1; k < stages; ++k) {
for (int j = 0; j < y0.length; ++j) {
T sum = yDotK[0][j].multiply(a[k - 1][0]);
for (int l = 1; l < k; ++l) {
sum = sum.add(yDotK[l][j].multiply(a[k - 1][l]));
}
yTmp[j] = y0[j].add(h.multiply(sum));
}
yDotK[k] = equations.computeDerivatives(t0.add(h.multiply(c[k - 1])), yTmp);
}
}
/** Apply internal weights of Butcher array, with corresponding times. Version with real Butcher array (non-Field).
* @param <T> the type of the field elements
* @param equations differential equations to integrate
* @param t0 initial time
* @param y0 initial value of the state vector at t0
* @param h step size
* @param a internal weights of Butcher array
* @param c times of Butcher array
* @param yDotK array where to store result
*/
static <T extends CalculusFieldElement<T>> void applyInternalButcherWeights(final FieldExpandableODE<T> equations,
final T t0, final T[] y0, final T h,
final double[][] a, final double[] c,
final T[][] yDotK) {
// create some internal working arrays
final int stages = c.length + 1;
final T[] yTmp = y0.clone();
for (int k = 1; k < stages; ++k) {
for (int j = 0; j < y0.length; ++j) {
T sum = yDotK[0][j].multiply(a[k - 1][0]);
for (int l = 1; l < k; ++l) {
sum = sum.add(yDotK[l][j].multiply(a[k - 1][l]));
}
yTmp[j] = y0[j].add(h.multiply(sum));
}
yDotK[k] = equations.computeDerivatives(t0.add(h.multiply(c[k - 1])), yTmp);
}
}
/** Apply external weights of Butcher array, assuming internal ones have been applied.
* @param <T> the type of the field elements
* @param yDotK output of stages
* @param y0 initial value of the state vector at t0
* @param h step size
* @param b external weights of Butcher array
* @return state vector
*/
static <T extends CalculusFieldElement<T>> T[] applyExternalButcherWeights(final T[] y0, final T[][] yDotK,
final T h, final T[] b) {
final T[] y = y0.clone();
final int stages = b.length;
for (int j = 0; j < y0.length; ++j) {
T sum = yDotK[0][j].multiply(b[0]);
for (int l = 1; l < stages; ++l) {
sum = sum.add(yDotK[l][j].multiply(b[l]));
}
y[j] = y[j].add(h.multiply(sum));
}
return y;
}
/** Apply external weights of Butcher array, assuming internal ones have been applied. Version with real Butcher
* array (non-Field version).
* @param <T> the type of the field elements
* @param yDotK output of stages
* @param y0 initial value of the state vector at t0
* @param h step size
* @param b external weights of Butcher array
* @return state vector
*/
static <T extends CalculusFieldElement<T>> T[] applyExternalButcherWeights(final T[] y0, final T[][] yDotK,
final T h, final double[] b) {
final T[] y = y0.clone();
final int stages = b.length;
for (int j = 0; j < y0.length; ++j) {
T sum = yDotK[0][j].multiply(b[0]);
for (int l = 1; l < stages; ++l) {
sum = sum.add(yDotK[l][j].multiply(b[l]));
}
y[j] = y[j].add(h.multiply(sum));
}
return y;
}
}