1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 package org.hipparchus.ode.nonstiff;
24
25
26 import org.hipparchus.CalculusFieldElement;
27 import org.hipparchus.Field;
28 import org.hipparchus.exception.MathIllegalArgumentException;
29 import org.hipparchus.exception.MathIllegalStateException;
30 import org.hipparchus.ode.AbstractFieldIntegrator;
31 import org.hipparchus.ode.FieldEquationsMapper;
32 import org.hipparchus.ode.FieldExpandableODE;
33 import org.hipparchus.ode.FieldODEState;
34 import org.hipparchus.ode.FieldODEStateAndDerivative;
35 import org.hipparchus.ode.nonstiff.interpolators.RungeKuttaFieldStateInterpolator;
36 import org.hipparchus.util.MathArrays;
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61 public abstract class FixedStepRungeKuttaFieldIntegrator<T extends CalculusFieldElement<T>>
62 extends AbstractFieldIntegrator<T> implements FieldExplicitRungeKuttaIntegrator<T> {
63
64
65 private final T[] c;
66
67
68 private final T[][] a;
69
70
71 private final T[] b;
72
73
74 private double[] realC = new double[0];
75
76
77 private double[][] realA = new double[0][];
78
79
80 private double[] realB = new double[0];
81
82
83 private final T step;
84
85
86 private boolean usingFieldCoefficients;
87
88
89
90
91
92
93
94
95 protected FixedStepRungeKuttaFieldIntegrator(final Field<T> field, final String name, final T step) {
96 super(field, name);
97 this.c = getC();
98 this.a = getA();
99 this.b = getB();
100 this.step = step.abs();
101 this.usingFieldCoefficients = false;
102 }
103
104
105
106
107 public T getDefaultStep() {
108 return this.step;
109 }
110
111
112
113
114
115
116 public void setUsingFieldCoefficients(boolean usingFieldCoefficients) {
117 this.usingFieldCoefficients = usingFieldCoefficients;
118 }
119
120
121 @Override
122 public boolean isUsingFieldCoefficients() {
123 return usingFieldCoefficients;
124 }
125
126
127 @Override
128 public int getNumberOfStages() {
129 return b.length;
130 }
131
132
133
134
135
136
137
138
139
140 protected abstract RungeKuttaFieldStateInterpolator<T> createInterpolator(boolean forward, T[][] yDotK,
141 FieldODEStateAndDerivative<T> globalPreviousState,
142 FieldODEStateAndDerivative<T> globalCurrentState,
143 FieldEquationsMapper<T> mapper);
144
145
146 @Override
147 protected FieldODEStateAndDerivative<T> initIntegration(FieldExpandableODE<T> eqn, FieldODEState<T> s0, T t) {
148 if (!isUsingFieldCoefficients()) {
149 realA = getRealA();
150 realB = getRealB();
151 realC = getRealC();
152 }
153 return super.initIntegration(eqn, s0, t);
154 }
155
156
157 @Override
158 public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
159 final FieldODEState<T> initialState, final T finalTime)
160 throws MathIllegalArgumentException, MathIllegalStateException {
161
162 sanityChecks(initialState, finalTime);
163 setStepStart(initIntegration(equations, initialState, finalTime));
164 final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0;
165
166
167 final int stages = getNumberOfStages();
168 final T[][] yDotK = MathArrays.buildArray(getField(), stages, -1);
169 MathArrays.buildArray(getField(), equations.getMapper().getTotalDimension());
170
171
172 if (forward) {
173 if (getStepStart().getTime().add(step).subtract(finalTime).getReal() >= 0) {
174 setStepSize(finalTime.subtract(getStepStart().getTime()));
175 } else {
176 setStepSize(step);
177 }
178 } else {
179 if (getStepStart().getTime().subtract(step).subtract(finalTime).getReal() <= 0) {
180 setStepSize(finalTime.subtract(getStepStart().getTime()));
181 } else {
182 setStepSize(step.negate());
183 }
184 }
185
186
187 setIsLastStep(false);
188 do {
189
190
191 final T[] y = getStepStart().getCompleteState();
192 yDotK[0] = getStepStart().getCompleteDerivative();
193
194
195 final T[] yTmp;
196 if (isUsingFieldCoefficients()) {
197 FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(), getStepStart().getTime(),
198 y, getStepSize(), a, c, yDotK);
199 yTmp = FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y, yDotK, getStepSize(), b);
200 } else {
201 FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(), getStepStart().getTime(),
202 y, getStepSize(), realA, realC, yDotK);
203 yTmp = FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y, yDotK, getStepSize(), realB);
204 }
205
206 incrementEvaluations(stages - 1);
207
208 final T stepEnd = getStepStart().getTime().add(getStepSize());
209 final T[] yDotTmp = computeDerivatives(stepEnd, yTmp);
210 final FieldODEStateAndDerivative<T> stateTmp = equations.getMapper().mapStateAndDerivative(stepEnd, yTmp, yDotTmp);
211
212
213 setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp, equations.getMapper()),
214 finalTime));
215
216 if (!isLastStep()) {
217
218
219 final T nextT = getStepStart().getTime().add(getStepSize());
220 final boolean nextIsLast = forward ?
221 (nextT.subtract(finalTime).getReal() >= 0) :
222 (nextT.subtract(finalTime).getReal() <= 0);
223 if (nextIsLast) {
224 setStepSize(finalTime.subtract(getStepStart().getTime()));
225 }
226 }
227
228 } while (!isLastStep());
229
230 final FieldODEStateAndDerivative<T> finalState = getStepStart();
231 setStepStart(null);
232 setStepSize(null);
233 return finalState;
234
235 }
236
237 }