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