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.nonstiff.interpolators.RungeKuttaStateInterpolator;
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 FixedStepRungeKuttaIntegrator extends AbstractIntegrator
55 implements ExplicitRungeKuttaIntegrator {
56
57
58 private final double[] c;
59
60
61 private final double[][] a;
62
63
64 private final double[] b;
65
66
67 private final double step;
68
69
70
71
72
73
74
75 protected FixedStepRungeKuttaIntegrator(final String name, final double step) {
76 super(name);
77 this.c = getC();
78 this.a = getA();
79 this.b = getB();
80 this.step = FastMath.abs(step);
81 }
82
83
84
85
86 public double getDefaultStep() {
87 return this.step;
88 }
89
90
91
92
93
94
95
96
97
98 protected abstract RungeKuttaStateInterpolator createInterpolator(boolean forward, double[][] yDotK,
99 ODEStateAndDerivative globalPreviousState,
100 ODEStateAndDerivative globalCurrentState,
101 EquationsMapper mapper);
102
103
104 @Override
105 public ODEStateAndDerivative integrate(final ExpandableODE equations,
106 final ODEState initialState, final double finalTime)
107 throws MathIllegalArgumentException, MathIllegalStateException {
108
109 sanityChecks(initialState, finalTime);
110 setStepStart(initIntegration(equations, initialState, finalTime));
111 final boolean forward = finalTime > initialState.getTime();
112
113
114 final int stages = c.length + 1;
115 double[] y = getStepStart().getCompleteState();
116 final double[][] yDotK = new double[stages][];
117 final double[] yTmp = new double[y.length];
118
119
120 if (forward) {
121 if (getStepStart().getTime() + step >= finalTime) {
122 setStepSize(finalTime - getStepStart().getTime());
123 } else {
124 setStepSize(step);
125 }
126 } else {
127 if (getStepStart().getTime() - step <= finalTime) {
128 setStepSize(finalTime - getStepStart().getTime());
129 } else {
130 setStepSize(-step);
131 }
132 }
133
134
135 setIsLastStep(false);
136 do {
137
138
139 y = getStepStart().getCompleteState();
140 yDotK[0] = getStepStart().getCompleteDerivative();
141
142
143 ExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(), getStepStart().getTime(), y,
144 getStepSize(), a, c, yDotK);
145
146 incrementEvaluations(stages - 1);
147
148
149 for (int j = 0; j < y.length; ++j) {
150 double sum = b[0] * yDotK[0][j];
151 for (int l = 1; l < stages; ++l) {
152 sum += b[l] * yDotK[l][j];
153 }
154 yTmp[j] = y[j] + getStepSize() * sum;
155 if (Double.isNaN(yTmp[j])) {
156 throw new MathIllegalStateException(LocalizedODEFormats.NAN_APPEARING_DURING_INTEGRATION,
157 getStepStart().getTime() + getStepSize());
158 }
159
160 }
161 final double stepEnd = getStepStart().getTime() + getStepSize();
162 final double[] yDotTmp = computeDerivatives(stepEnd, yTmp);
163 final ODEStateAndDerivative stateTmp =
164 equations.getMapper().mapStateAndDerivative(stepEnd, yTmp, yDotTmp);
165
166
167 System.arraycopy(yTmp, 0, y, 0, y.length);
168 setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp,
169 equations.getMapper()),
170 finalTime));
171
172 if (!isLastStep()) {
173
174
175 final double nextT = getStepStart().getTime() + getStepSize();
176 final boolean nextIsLast = forward ? (nextT >= finalTime) : (nextT <= finalTime);
177 if (nextIsLast) {
178 setStepSize(finalTime - getStepStart().getTime());
179 }
180 }
181
182 } while (!isLastStep());
183
184 final ODEStateAndDerivative finalState = getStepStart();
185 setStepStart(null);
186 setStepSize(Double.NaN);
187 return finalState;
188
189 }
190
191 }