1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math.ode.nonstiff;
19
20 import org.apache.commons.math.ode.AbstractIntegrator;
21 import org.apache.commons.math.ode.DerivativeException;
22 import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
23 import org.apache.commons.math.ode.IntegratorException;
24
25
26
27
28
29
30
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 AdaptiveStepsizeIntegrator
55 extends AbstractIntegrator {
56
57
58
59
60
61
62
63
64
65
66
67
68 public AdaptiveStepsizeIntegrator(final String name,
69 final double minStep, final double maxStep,
70 final double scalAbsoluteTolerance,
71 final double scalRelativeTolerance) {
72
73 super(name);
74
75 this.minStep = Math.abs(minStep);
76 this.maxStep = Math.abs(maxStep);
77 this.initialStep = -1.0;
78
79 this.scalAbsoluteTolerance = scalAbsoluteTolerance;
80 this.scalRelativeTolerance = scalRelativeTolerance;
81 this.vecAbsoluteTolerance = null;
82 this.vecRelativeTolerance = null;
83
84 resetInternalState();
85
86 }
87
88
89
90
91
92
93
94
95
96
97
98 public AdaptiveStepsizeIntegrator(final String name,
99 final double minStep, final double maxStep,
100 final double[] vecAbsoluteTolerance,
101 final double[] vecRelativeTolerance) {
102
103 super(name);
104
105 this.minStep = minStep;
106 this.maxStep = maxStep;
107 this.initialStep = -1.0;
108
109 this.scalAbsoluteTolerance = 0;
110 this.scalRelativeTolerance = 0;
111 this.vecAbsoluteTolerance = vecAbsoluteTolerance.clone();
112 this.vecRelativeTolerance = vecRelativeTolerance.clone();
113
114 resetInternalState();
115
116 }
117
118
119
120
121
122
123
124
125
126
127
128
129 public void setInitialStepSize(final double initialStepSize) {
130 if ((initialStepSize < minStep) || (initialStepSize > maxStep)) {
131 initialStep = -1.0;
132 } else {
133 initialStep = initialStepSize;
134 }
135 }
136
137
138
139
140
141
142
143
144
145 @Override
146 protected void sanityChecks(final FirstOrderDifferentialEquations equations,
147 final double t0, final double[] y0,
148 final double t, final double[] y)
149 throws IntegratorException {
150
151 super.sanityChecks(equations, t0, y0, t, y);
152
153 if ((vecAbsoluteTolerance != null) && (vecAbsoluteTolerance.length != y0.length)) {
154 throw new IntegratorException(
155 "dimensions mismatch: state vector has dimension {0}," +
156 " absolute tolerance vector has dimension {1}",
157 y0.length, vecAbsoluteTolerance.length);
158 }
159
160 if ((vecRelativeTolerance != null) && (vecRelativeTolerance.length != y0.length)) {
161 throw new IntegratorException(
162 "dimensions mismatch: state vector has dimension {0}," +
163 " relative tolerance vector has dimension {1}",
164 y0.length, vecRelativeTolerance.length);
165 }
166
167 }
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183 public double initializeStep(final FirstOrderDifferentialEquations equations,
184 final boolean forward, final int order, final double[] scale,
185 final double t0, final double[] y0, final double[] yDot0,
186 final double[] y1, final double[] yDot1)
187 throws DerivativeException {
188
189 if (initialStep > 0) {
190
191 return forward ? initialStep : -initialStep;
192 }
193
194
195
196 double ratio;
197 double yOnScale2 = 0;
198 double yDotOnScale2 = 0;
199 for (int j = 0; j < y0.length; ++j) {
200 ratio = y0[j] / scale[j];
201 yOnScale2 += ratio * ratio;
202 ratio = yDot0[j] / scale[j];
203 yDotOnScale2 += ratio * ratio;
204 }
205
206 double h = ((yOnScale2 < 1.0e-10) || (yDotOnScale2 < 1.0e-10)) ?
207 1.0e-6 : (0.01 * Math.sqrt(yOnScale2 / yDotOnScale2));
208 if (! forward) {
209 h = -h;
210 }
211
212
213 for (int j = 0; j < y0.length; ++j) {
214 y1[j] = y0[j] + h * yDot0[j];
215 }
216 computeDerivatives(t0 + h, y1, yDot1);
217
218
219 double yDDotOnScale = 0;
220 for (int j = 0; j < y0.length; ++j) {
221 ratio = (yDot1[j] - yDot0[j]) / scale[j];
222 yDDotOnScale += ratio * ratio;
223 }
224 yDDotOnScale = Math.sqrt(yDDotOnScale) / h;
225
226
227
228 final double maxInv2 = Math.max(Math.sqrt(yDotOnScale2), yDDotOnScale);
229 final double h1 = (maxInv2 < 1.0e-15) ?
230 Math.max(1.0e-6, 0.001 * Math.abs(h)) :
231 Math.pow(0.01 / maxInv2, 1.0 / order);
232 h = Math.min(100.0 * Math.abs(h), h1);
233 h = Math.max(h, 1.0e-12 * Math.abs(t0));
234 if (h < getMinStep()) {
235 h = getMinStep();
236 }
237 if (h > getMaxStep()) {
238 h = getMaxStep();
239 }
240 if (! forward) {
241 h = -h;
242 }
243
244 return h;
245
246 }
247
248
249
250
251
252
253
254
255
256
257 protected double filterStep(final double h, final boolean forward, final boolean acceptSmall)
258 throws IntegratorException {
259
260 double filteredH = h;
261 if (Math.abs(h) < minStep) {
262 if (acceptSmall) {
263 filteredH = forward ? minStep : -minStep;
264 } else {
265 throw new IntegratorException(
266 "minimal step size ({0}) reached, integration needs {1}",
267 minStep, Math.abs(h));
268 }
269 }
270
271 if (filteredH > maxStep) {
272 filteredH = maxStep;
273 } else if (filteredH < -maxStep) {
274 filteredH = -maxStep;
275 }
276
277 return filteredH;
278
279 }
280
281
282 public abstract double integrate (FirstOrderDifferentialEquations equations,
283 double t0, double[] y0,
284 double t, double[] y)
285 throws DerivativeException, IntegratorException;
286
287
288 @Override
289 public double getCurrentStepStart() {
290 return stepStart;
291 }
292
293
294 protected void resetInternalState() {
295 stepStart = Double.NaN;
296 stepSize = Math.sqrt(minStep * maxStep);
297 }
298
299
300
301
302 public double getMinStep() {
303 return minStep;
304 }
305
306
307
308
309 public double getMaxStep() {
310 return maxStep;
311 }
312
313
314 private double minStep;
315
316
317 private double maxStep;
318
319
320 private double initialStep;
321
322
323 protected double scalAbsoluteTolerance;
324
325
326 protected double scalRelativeTolerance;
327
328
329 protected double[] vecAbsoluteTolerance;
330
331
332 protected double[] vecRelativeTolerance;
333
334 }