1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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.ode.OrdinaryDifferentialEquation;
30 import org.hipparchus.util.FastMath;
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54 public abstract class RungeKuttaIntegrator extends AbstractIntegrator implements ButcherArrayProvider {
55
56
57 private final double[] c;
58
59
60 private final double[][] a;
61
62
63 private final double[] b;
64
65
66 private final double step;
67
68
69
70
71
72
73
74 protected RungeKuttaIntegrator(final String name, final double step) {
75 super(name);
76 this.c = getC();
77 this.a = getA();
78 this.b = getB();
79 this.step = FastMath.abs(step);
80 }
81
82
83
84
85 public double getDefaultStep() {
86 return this.step;
87 }
88
89
90
91
92
93
94
95
96
97 protected abstract RungeKuttaStateInterpolator createInterpolator(boolean forward, double[][] yDotK,
98 ODEStateAndDerivative globalPreviousState,
99 ODEStateAndDerivative globalCurrentState,
100 EquationsMapper mapper);
101
102
103 @Override
104 public ODEStateAndDerivative integrate(final ExpandableODE equations,
105 final ODEState initialState, final double finalTime)
106 throws MathIllegalArgumentException, MathIllegalStateException {
107
108 sanityChecks(initialState, finalTime);
109 setStepStart(initIntegration(equations, initialState, finalTime));
110 final boolean forward = finalTime > initialState.getTime();
111
112
113 final int stages = c.length + 1;
114 double[] y = getStepStart().getCompleteState();
115 final double[][] yDotK = new double[stages][];
116 final double[] yTmp = new double[y.length];
117
118
119 if (forward) {
120 if (getStepStart().getTime() + step >= finalTime) {
121 setStepSize(finalTime - getStepStart().getTime());
122 } else {
123 setStepSize(step);
124 }
125 } else {
126 if (getStepStart().getTime() - step <= finalTime) {
127 setStepSize(finalTime - getStepStart().getTime());
128 } else {
129 setStepSize(-step);
130 }
131 }
132
133
134 setIsLastStep(false);
135 do {
136
137
138 y = getStepStart().getCompleteState();
139 yDotK[0] = getStepStart().getCompleteDerivative();
140
141
142 for (int k = 1; k < stages; ++k) {
143
144 for (int j = 0; j < y.length; ++j) {
145 double sum = a[k-1][0] * yDotK[0][j];
146 for (int l = 1; l < k; ++l) {
147 sum += a[k-1][l] * yDotK[l][j];
148 }
149 yTmp[j] = y[j] + getStepSize() * sum;
150 }
151
152 yDotK[k] = computeDerivatives(getStepStart().getTime() + c[k-1] * getStepSize(), yTmp);
153
154 }
155
156
157 for (int j = 0; j < y.length; ++j) {
158 double sum = b[0] * yDotK[0][j];
159 for (int l = 1; l < stages; ++l) {
160 sum += b[l] * yDotK[l][j];
161 }
162 yTmp[j] = y[j] + getStepSize() * sum;
163 if (Double.isNaN(yTmp[j])) {
164 throw new MathIllegalStateException(LocalizedODEFormats.NAN_APPEARING_DURING_INTEGRATION,
165 getStepStart().getTime() + getStepSize());
166 }
167
168 }
169 final double stepEnd = getStepStart().getTime() + getStepSize();
170 final double[] yDotTmp = computeDerivatives(stepEnd, yTmp);
171 final ODEStateAndDerivative stateTmp =
172 equations.getMapper().mapStateAndDerivative(stepEnd, yTmp, yDotTmp);
173
174
175 System.arraycopy(yTmp, 0, y, 0, y.length);
176 setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp,
177 equations.getMapper()),
178 finalTime));
179
180 if (!isLastStep()) {
181
182
183 final double nextT = getStepStart().getTime() + getStepSize();
184 final boolean nextIsLast = forward ? (nextT >= finalTime) : (nextT <= finalTime);
185 if (nextIsLast) {
186 setStepSize(finalTime - getStepStart().getTime());
187 }
188 }
189
190 } while (!isLastStep());
191
192 final ODEStateAndDerivative finalState = getStepStart();
193 setStepStart(null);
194 setStepSize(Double.NaN);
195 return finalState;
196
197 }
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224 public double[] singleStep(final OrdinaryDifferentialEquation equations,
225 final double t0, final double[] y0, final double t) {
226
227
228 final double[] y = y0.clone();
229 final int stages = c.length + 1;
230 final double[][] yDotK = new double[stages][];
231 final double[] yTmp = y0.clone();
232
233
234 final double h = t - t0;
235 yDotK[0] = equations.computeDerivatives(t0, y);
236
237
238 for (int k = 1; k < stages; ++k) {
239
240 for (int j = 0; j < y0.length; ++j) {
241 double sum = a[k-1][0] * yDotK[0][j];
242 for (int l = 1; l < k; ++l) {
243 sum += a[k-1][l] * yDotK[l][j];
244 }
245 yTmp[j] = y[j] + h * sum;
246 }
247
248 yDotK[k] = equations.computeDerivatives(t0 + c[k-1] * h, yTmp);
249
250 }
251
252
253 for (int j = 0; j < y0.length; ++j) {
254 double sum = b[0] * yDotK[0][j];
255 for (int l = 1; l < stages; ++l) {
256 sum += b[l] * yDotK[l][j];
257 }
258 y[j] += h * sum;
259 }
260
261 return y;
262
263 }
264
265 }