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.exception.MathIllegalArgumentException;
22 import org.hipparchus.exception.MathIllegalStateException;
23 import org.hipparchus.ode.AbstractIntegrator;
24 import org.hipparchus.ode.EquationsMapper;
25 import org.hipparchus.ode.ExpandableODE;
26 import org.hipparchus.ode.LocalizedODEFormats;
27 import org.hipparchus.ode.ODEState;
28 import org.hipparchus.ode.ODEStateAndDerivative;
29 import org.hipparchus.util.FastMath;
30
31 /**
32 * This class implements the common part of all fixed step Runge-Kutta
33 * integrators for Ordinary Differential Equations.
34 *
35 * <p>These methods are explicit Runge-Kutta methods, their Butcher
36 * arrays are as follows :</p>
37 * <pre>
38 * 0 |
39 * c2 | a21
40 * c3 | a31 a32
41 * ... | ...
42 * cs | as1 as2 ... ass-1
43 * |--------------------------
44 * | b1 b2 ... bs-1 bs
45 * </pre>
46 *
47 * @see EulerIntegrator
48 * @see ClassicalRungeKuttaIntegrator
49 * @see GillIntegrator
50 * @see MidpointIntegrator
51 */
52
53 public abstract class RungeKuttaIntegrator extends AbstractIntegrator implements ExplicitRungeKuttaIntegrator {
54
55 /** Time steps from Butcher array (without the first zero). */
56 private final double[] c;
57
58 /** Internal weights from Butcher array (without the first empty row). */
59 private final double[][] a;
60
61 /** External weights for the high order method from Butcher array. */
62 private final double[] b;
63
64 /** Integration step. */
65 private final double step;
66
67 /** Simple constructor.
68 * Build a Runge-Kutta integrator with the given
69 * step. The default step handler does nothing.
70 * @param name name of the method
71 * @param step integration step
72 */
73 protected RungeKuttaIntegrator(final String name, final double step) {
74 super(name);
75 this.c = getC();
76 this.a = getA();
77 this.b = getB();
78 this.step = FastMath.abs(step);
79 }
80
81 /** Getter for the default, positive step-size assigned at constructor level.
82 * @return step
83 */
84 public double getDefaultStep() {
85 return this.step;
86 }
87
88 /** Create an interpolator.
89 * @param forward integration direction indicator
90 * @param yDotK slopes at the intermediate points
91 * @param globalPreviousState start of the global step
92 * @param globalCurrentState end of the global step
93 * @param mapper equations mapper for the all equations
94 * @return external weights for the high order method from Butcher array
95 */
96 protected abstract RungeKuttaStateInterpolator createInterpolator(boolean forward, double[][] yDotK,
97 ODEStateAndDerivative globalPreviousState,
98 ODEStateAndDerivative globalCurrentState,
99 EquationsMapper mapper);
100
101 /** {@inheritDoc} */
102 @Override
103 public ODEStateAndDerivative integrate(final ExpandableODE equations,
104 final ODEState initialState, final double finalTime)
105 throws MathIllegalArgumentException, MathIllegalStateException {
106
107 sanityChecks(initialState, finalTime);
108 setStepStart(initIntegration(equations, initialState, finalTime));
109 final boolean forward = finalTime > initialState.getTime();
110
111 // create some internal working arrays
112 final int stages = c.length + 1;
113 double[] y = getStepStart().getCompleteState();
114 final double[][] yDotK = new double[stages][];
115 final double[] yTmp = new double[y.length];
116
117 // set up integration control objects
118 if (forward) {
119 if (getStepStart().getTime() + step >= finalTime) {
120 setStepSize(finalTime - getStepStart().getTime());
121 } else {
122 setStepSize(step);
123 }
124 } else {
125 if (getStepStart().getTime() - step <= finalTime) {
126 setStepSize(finalTime - getStepStart().getTime());
127 } else {
128 setStepSize(-step);
129 }
130 }
131
132 // main integration loop
133 setIsLastStep(false);
134 do {
135
136 // first stage
137 y = getStepStart().getCompleteState();
138 yDotK[0] = getStepStart().getCompleteDerivative();
139
140 // next stages
141 ExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(), getStepStart().getTime(), y,
142 getStepSize(), a, c, yDotK);
143
144 incrementEvaluations(stages - 1);
145
146 // estimate the state at the end of the step
147 for (int j = 0; j < y.length; ++j) {
148 double sum = b[0] * yDotK[0][j];
149 for (int l = 1; l < stages; ++l) {
150 sum += b[l] * yDotK[l][j];
151 }
152 yTmp[j] = y[j] + getStepSize() * sum;
153 if (Double.isNaN(yTmp[j])) {
154 throw new MathIllegalStateException(LocalizedODEFormats.NAN_APPEARING_DURING_INTEGRATION,
155 getStepStart().getTime() + getStepSize());
156 }
157
158 }
159 final double stepEnd = getStepStart().getTime() + getStepSize();
160 final double[] yDotTmp = computeDerivatives(stepEnd, yTmp);
161 final ODEStateAndDerivative stateTmp =
162 equations.getMapper().mapStateAndDerivative(stepEnd, yTmp, yDotTmp);
163
164 // discrete events handling
165 System.arraycopy(yTmp, 0, y, 0, y.length);
166 setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp,
167 equations.getMapper()),
168 finalTime));
169
170 if (!isLastStep()) {
171
172 // stepsize control for next step
173 final double nextT = getStepStart().getTime() + getStepSize();
174 final boolean nextIsLast = forward ? (nextT >= finalTime) : (nextT <= finalTime);
175 if (nextIsLast) {
176 setStepSize(finalTime - getStepStart().getTime());
177 }
178 }
179
180 } while (!isLastStep());
181
182 final ODEStateAndDerivative finalState = getStepStart();
183 setStepStart(null);
184 setStepSize(Double.NaN);
185 return finalState;
186
187 }
188
189 }