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 package org.hipparchus.ode.nonstiff;
18
19 import org.hipparchus.ode.ExpandableODE;
20 import org.hipparchus.ode.ODEIntegrator;
21 import org.hipparchus.ode.ODEState;
22 import org.hipparchus.ode.OrdinaryDifferentialEquation;
23
24
25 /**
26 * This interface implements the part of Runge-Kutta
27 * integrators for Ordinary Differential Equations
28 * common to fixed- and adaptive steps.
29 *
30 * <p>These methods are explicit Runge-Kutta methods, their Butcher
31 * arrays are as follows :</p>
32 * <pre>
33 * 0 |
34 * c2 | a21
35 * c3 | a31 a32
36 * ... | ...
37 * cs | as1 as2 ... ass-1
38 * |--------------------------
39 * | b1 b2 ... bs-1 bs
40 * </pre>
41 *
42 * @see ButcherArrayProvider
43 * @see RungeKuttaIntegrator
44 * @see EmbeddedRungeKuttaIntegrator
45 * @since 3.1
46 */
47
48 public interface ExplicitRungeKuttaIntegrator extends ButcherArrayProvider, ODEIntegrator {
49
50 /**
51 * Getter for the number of stages corresponding to the Butcher array.
52 *
53 * @return number of stages
54 */
55 default int getNumberOfStages() {
56 return getB().length;
57 }
58
59 /** Fast computation of a single step of ODE integration.
60 * <p>This method is intended for the limited use case of
61 * very fast computation of only one step without using any of the
62 * rich features of general integrators that may take some time
63 * to set up (i.e. no step handlers, no events handlers, no additional
64 * states, no interpolators, no error control, no evaluations count,
65 * no sanity checks ...). It handles the strict minimum of computation,
66 * so it can be embedded in outer loops.</p>
67 * <p>
68 * This method is <em>not</em> used at all by the {@link #integrate(ExpandableODE, ODEState, double)}
69 * method. It also completely ignores the step set at construction time, and
70 * uses only a single step to go from {@code t0} to {@code t}.
71 * </p>
72 * <p>
73 * As this method does not use any of the state-dependent features of the integrator,
74 * it should be reasonably thread-safe <em>if and only if</em> the provided differential
75 * equations are themselves thread-safe.
76 * </p>
77 * @param equations differential equations to integrate
78 * @param t0 initial time
79 * @param y0 initial value of the state vector at t0
80 * @param t target time for the integration
81 * (can be set to a value smaller than {@code t0} for backward integration)
82 * @return state vector at {@code t}
83 */
84 default double[] singleStep(final OrdinaryDifferentialEquation equations, final double t0, final double[] y0,
85 final double t) {
86
87 // create some internal working arrays
88 final int stages = getNumberOfStages();
89 final double[][] yDotK = new double[stages][];
90
91 // first stage
92 final double h = t - t0;
93 final ExpandableODE expandableODE = new ExpandableODE(equations);
94 yDotK[0] = expandableODE.computeDerivatives(t0, y0);
95
96 // next stages
97 ExplicitRungeKuttaIntegrator.applyInternalButcherWeights(expandableODE, t0, y0, h, getA(), getC(), yDotK);
98
99 // estimate the state at the end of the step
100 return ExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y0, yDotK, h, getB());
101
102 }
103
104 /**
105 * Apply internal weights of Butcher array, with corresponding times.
106 * @param equations differential equations to integrate
107 * @param t0 initial time
108 * @param y0 initial value of the state vector at t0
109 * @param h step size
110 * @param a internal weights of Butcher array
111 * @param c times of Butcher array
112 * @param yDotK array where to store result
113 */
114 static void applyInternalButcherWeights(final ExpandableODE equations, final double t0, final double[] y0,
115 final double h, final double[][] a, final double[] c,
116 final double[][] yDotK) {
117 // create some internal working arrays
118 final int stages = c.length + 1;
119 final double[] yTmp = y0.clone();
120
121 for (int k = 1; k < stages; ++k) {
122
123 for (int j = 0; j < y0.length; ++j) {
124 double sum = yDotK[0][j] * a[k - 1][0];
125 for (int l = 1; l < k; ++l) {
126 sum += yDotK[l][j] * a[k - 1][l];
127 }
128 yTmp[j] = y0[j] + h * sum;
129 }
130
131 yDotK[k] = equations.computeDerivatives(t0 + h * c[k - 1], yTmp);
132 }
133 }
134
135 /** Apply external weights of Butcher array, assuming internal ones have been applied.
136 * @param yDotK output of stages
137 * @param y0 initial value of the state vector at t0
138 * @param h step size
139 * @param b external weights of Butcher array
140 * @return state vector
141 */
142 static double[] applyExternalButcherWeights(final double[] y0, final double[][] yDotK, final double h,
143 final double[] b) {
144 final double[] y = y0.clone();
145 final int stages = b.length;
146 for (int j = 0; j < y0.length; ++j) {
147 double sum = yDotK[0][j] * b[0];
148 for (int l = 1; l < stages; ++l) {
149 sum += yDotK[l][j] * b[l];
150 }
151 y[j] += h * sum;
152 }
153 return y;
154 }
155 }