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 import org.hipparchus.CalculusFieldElement;
21 import org.hipparchus.Field;
22 import org.hipparchus.analysis.differentiation.DSFactory;
23 import org.hipparchus.analysis.differentiation.DerivativeStructure;
24 import org.hipparchus.analysis.differentiation.Gradient;
25 import org.hipparchus.analysis.solvers.BracketedRealFieldUnivariateSolver;
26 import org.hipparchus.analysis.solvers.FieldBracketingNthOrderBrentSolver;
27 import org.hipparchus.complex.Complex;
28 import org.hipparchus.complex.ComplexField;
29 import org.hipparchus.exception.LocalizedCoreFormats;
30 import org.hipparchus.exception.MathIllegalArgumentException;
31 import org.hipparchus.exception.MathIllegalStateException;
32 import org.hipparchus.ode.FieldExpandableODE;
33 import org.hipparchus.ode.FieldODEIntegrator;
34 import org.hipparchus.ode.FieldODEState;
35 import org.hipparchus.ode.FieldODEStateAndDerivative;
36 import org.hipparchus.ode.FieldOrdinaryDifferentialEquation;
37 import org.hipparchus.ode.FieldSecondaryODE;
38 import org.hipparchus.ode.LocalizedODEFormats;
39 import org.hipparchus.ode.TestFieldProblem1;
40 import org.hipparchus.ode.TestFieldProblem2;
41 import org.hipparchus.ode.TestFieldProblem3;
42 import org.hipparchus.ode.TestFieldProblem4;
43 import org.hipparchus.ode.TestFieldProblem5;
44 import org.hipparchus.ode.TestFieldProblem6;
45 import org.hipparchus.ode.TestFieldProblemAbstract;
46 import org.hipparchus.ode.TestFieldProblemHandler;
47 import org.hipparchus.ode.events.Action;
48 import org.hipparchus.ode.events.FieldAdaptableInterval;
49 import org.hipparchus.ode.events.FieldODEEventDetector;
50 import org.hipparchus.ode.events.FieldODEEventHandler;
51 import org.hipparchus.ode.sampling.FieldODEStateInterpolator;
52 import org.hipparchus.ode.sampling.FieldODEStepHandler;
53 import org.hipparchus.ode.sampling.StepInterpolatorTestUtils;
54 import org.hipparchus.util.FastMath;
55 import org.hipparchus.util.MathArrays;
56 import org.junit.jupiter.api.Test;
57
58 import java.lang.reflect.Array;
59 import java.lang.reflect.Constructor;
60 import java.lang.reflect.InvocationTargetException;
61
62 import static org.junit.jupiter.api.Assertions.assertEquals;
63 import static org.junit.jupiter.api.Assertions.assertNotNull;
64 import static org.junit.jupiter.api.Assertions.assertSame;
65 import static org.junit.jupiter.api.Assertions.assertTrue;
66 import static org.junit.jupiter.api.Assertions.fail;
67
68 public abstract class RungeKuttaFieldIntegratorAbstractTest {
69
70 protected abstract <T extends CalculusFieldElement<T>> FixedStepRungeKuttaFieldIntegrator<T>
71 createIntegrator(Field<T> field, T step);
72
73 @Test
74 public abstract void testNonFieldIntegratorConsistency();
75
76 protected <T extends CalculusFieldElement<T>> void doTestNonFieldIntegratorConsistency(final Field<T> field) {
77 try {
78
79
80 FixedStepRungeKuttaFieldIntegrator<T> fieldIntegrator = createIntegrator(field, field.getZero().add(1));
81 T[][] fieldA = fieldIntegrator.getA();
82 T[] fieldB = fieldIntegrator.getB();
83 T[] fieldC = fieldIntegrator.getC();
84
85 String fieldName = fieldIntegrator.getClass().getName();
86 String regularName = fieldName.replaceAll("Field", "");
87
88
89 @SuppressWarnings("unchecked")
90 Constructor<FixedStepRungeKuttaIntegrator> constructor =
91 (Constructor<FixedStepRungeKuttaIntegrator>) Class.forName(regularName).getConstructor(Double.TYPE);
92 final FixedStepRungeKuttaIntegrator regularIntegrator =
93 constructor.newInstance(1.0);
94 double[][] regularA = regularIntegrator.getA();
95 double[] regularB = regularIntegrator.getB();
96 double[] regularC = regularIntegrator.getC();
97
98 assertEquals(regularA.length, fieldA.length);
99 for (int i = 0; i < regularA.length; ++i) {
100 checkArray(regularA[i], fieldA[i]);
101 }
102 checkArray(regularB, fieldB);
103 checkArray(regularC, fieldC);
104
105 } catch (ClassNotFoundException | IllegalAccessException | IllegalArgumentException |
106 SecurityException | NoSuchMethodException | InvocationTargetException |
107 InstantiationException e) {
108 fail(e.getLocalizedMessage());
109 }
110 }
111
112 private <T extends CalculusFieldElement<T>> void checkArray(double[] regularArray, T[] fieldArray) {
113 assertEquals(regularArray.length, fieldArray.length);
114 for (int i = 0; i < regularArray.length; ++i) {
115 if (regularArray[i] == 0) {
116 assertEquals(0.0, fieldArray[i].getReal());
117 } else {
118 assertEquals(regularArray[i], fieldArray[i].getReal(), FastMath.ulp(regularArray[i]));
119 }
120 }
121 }
122
123 @Test
124 public abstract void testMissedEndEvent();
125
126 protected <T extends CalculusFieldElement<T>> void doTestMissedEndEvent(final Field<T> field,
127 final double epsilonT, final double epsilonY)
128 throws MathIllegalArgumentException, MathIllegalStateException {
129 final T t0 = field.getZero().add(1878250320.0000029);
130 final T tEvent = field.getZero().add(1878250379.9999986);
131 final T[] k = MathArrays.buildArray(field, 3);
132 k[0] = field.getZero().add(1.0e-4);
133 k[1] = field.getZero().add(1.0e-5);
134 k[2] = field.getZero().add(1.0e-6);
135 FieldOrdinaryDifferentialEquation<T> ode = new FieldOrdinaryDifferentialEquation<T>() {
136
137 public int getDimension() {
138 return k.length;
139 }
140
141 public T[] computeDerivatives(T t, T[] y) {
142 T[] yDot = MathArrays.buildArray(field, k.length);
143 for (int i = 0; i < y.length; ++i) {
144 yDot[i] = k[i].multiply(y[i]);
145 }
146 return yDot;
147 }
148 };
149
150 FixedStepRungeKuttaFieldIntegrator<T> integrator = createIntegrator(field, field.getZero().add(60.0));
151
152 T[] y0 = MathArrays.buildArray(field, k.length);
153 for (int i = 0; i < y0.length; ++i) {
154 y0[i] = field.getOne().add(i);
155 }
156
157 FieldODEStateAndDerivative<T> result = integrator.integrate(new FieldExpandableODE<T>(ode),
158 new FieldODEState<T>(t0, y0),
159 tEvent);
160 assertEquals(tEvent.getReal(), result.getTime().getReal(), epsilonT);
161 T[] y = result.getPrimaryState();
162 for (int i = 0; i < y.length; ++i) {
163 assertEquals(y0[i].multiply(k[i].multiply(result.getTime().subtract(t0)).exp()).getReal(),
164 y[i].getReal(),
165 epsilonY);
166 }
167
168 integrator.addEventDetector(new FieldODEEventDetector<T>() {
169 public FieldAdaptableInterval<T> getMaxCheckInterval() {
170 return (s, isForward) -> Double.POSITIVE_INFINITY;
171 }
172 public int getMaxIterationCount() {
173 return 100;
174 }
175 public BracketedRealFieldUnivariateSolver<T> getSolver() {
176 return new FieldBracketingNthOrderBrentSolver<T>(field.getZero(),
177 field.getZero().newInstance(1.0e-20),
178 field.getZero(),
179 5);
180 }
181 public T g(FieldODEStateAndDerivative<T> state) {
182 return state.getTime().subtract(tEvent);
183 }
184 public FieldODEEventHandler<T> getHandler() {
185 return (state, detector, increasing) -> {
186 assertEquals(tEvent.getReal(), state.getTime().getReal(), epsilonT);
187 return Action.CONTINUE;
188 };
189 }
190 });
191 result = integrator.integrate(new FieldExpandableODE<T>(ode),
192 new FieldODEState<T>(t0, y0),
193 tEvent.add(120));
194 assertEquals(tEvent.add(120).getReal(), result.getTime().getReal(), epsilonT);
195 y = result.getPrimaryState();
196 for (int i = 0; i < y.length; ++i) {
197 assertEquals(y0[i].multiply(k[i].multiply(result.getTime().subtract(t0)).exp()).getReal(),
198 y[i].getReal(),
199 epsilonY);
200 }
201
202 }
203
204 @Test
205 public abstract void testSanityChecks();
206
207 protected <T extends CalculusFieldElement<T>> void doTestSanityChecks(Field<T> field)
208 throws MathIllegalArgumentException, MathIllegalStateException {
209 FixedStepRungeKuttaFieldIntegrator<T> integrator = createIntegrator(field, field.getZero().add(0.01));
210 try {
211 TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
212 integrator.integrate(new FieldExpandableODE<T>(pb),
213 new FieldODEState<T>(field.getZero(), MathArrays.buildArray(field, pb.getDimension() + 10)),
214 field.getOne());
215 fail("an exception should have been thrown");
216 } catch(MathIllegalArgumentException ie) {
217 assertEquals(LocalizedCoreFormats.DIMENSIONS_MISMATCH, ie.getSpecifier());
218 }
219 try {
220 TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
221 integrator.integrate(new FieldExpandableODE<T>(pb),
222 new FieldODEState<T>(field.getZero(), MathArrays.buildArray(field, pb.getDimension())),
223 field.getZero());
224 fail("an exception should have been thrown");
225 } catch(MathIllegalArgumentException ie) {
226 assertEquals(LocalizedODEFormats.TOO_SMALL_INTEGRATION_INTERVAL, ie.getSpecifier());
227 }
228 }
229
230 @Test
231 public abstract void testDecreasingSteps();
232
233 protected <T extends CalculusFieldElement<T>> void doTestDecreasingSteps(Field<T> field,
234 final double safetyValueFactor,
235 final double safetyTimeFactor,
236 final double epsilonT)
237 throws MathIllegalArgumentException, MathIllegalStateException {
238
239 @SuppressWarnings("unchecked")
240 TestFieldProblemAbstract<T>[] allProblems =
241 (TestFieldProblemAbstract<T>[]) Array.newInstance(TestFieldProblemAbstract.class, 6);
242 allProblems[0] = new TestFieldProblem1<T>(field);
243 allProblems[1] = new TestFieldProblem2<T>(field);
244 allProblems[2] = new TestFieldProblem3<T>(field);
245 allProblems[3] = new TestFieldProblem4<T>(field);
246 allProblems[4] = new TestFieldProblem5<T>(field);
247 allProblems[5] = new TestFieldProblem6<T>(field);
248 for (TestFieldProblemAbstract<T> pb : allProblems) {
249
250 T previousValueError = null;
251 T previousTimeError = null;
252 for (int i = 4; i < 10; ++i) {
253
254 T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(FastMath.pow(2.0, -i));
255
256 FixedStepRungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step);
257 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
258 integ.addStepHandler(handler);
259 final double maxCheck = Double.POSITIVE_INFINITY;
260 final T eventTol = step.multiply(1.0e-6);
261 FieldODEEventDetector<T>[] functions = pb.getEventDetectors(maxCheck, eventTol, 1000);
262 for (int l = 0; l < functions.length; ++l) {
263 integ.addEventDetector(functions[l]);
264 }
265 assertEquals(functions.length, integ.getEventDetectors().size());
266 FieldODEStateAndDerivative<T> stop = integ.integrate(new FieldExpandableODE<T>(pb),
267 pb.getInitialState(),
268 pb.getFinalTime());
269 if (functions.length == 0) {
270 assertEquals(pb.getFinalTime().getReal(), stop.getTime().getReal(), epsilonT);
271 }
272
273 T error = handler.getMaximalValueError();
274 if (i > 4) {
275 assertTrue(error.subtract(previousValueError.abs().multiply(safetyValueFactor)).getReal() < 0);
276 }
277 previousValueError = error;
278
279 T timeError = handler.getMaximalTimeError();
280 if (i > 4) {
281
282 T timeTol = max(eventTol, previousTimeError.abs().multiply(safetyTimeFactor));
283 assertTrue(timeError.subtract(timeTol).getReal() <= 0);
284 }
285 previousTimeError = timeError;
286
287 integ.clearEventDetectors();
288 assertEquals(0, integ.getEventDetectors().size());
289 }
290
291 }
292
293 }
294
295
296
297
298
299
300
301
302 private <T extends CalculusFieldElement<T>> T max(T a, T b) {
303 return a.getReal() > b.getReal() ? a : b;
304 }
305
306 @Test
307 public abstract void testSmallStep();
308
309 protected <T extends CalculusFieldElement<T>> void doTestSmallStep(Field<T> field,
310 final double epsilonLast,
311 final double epsilonMaxValue,
312 final double epsilonMaxTime,
313 final String name)
314 throws MathIllegalArgumentException, MathIllegalStateException {
315
316 TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
317 T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.001);
318
319 FixedStepRungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step);
320 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
321 integ.addStepHandler(handler);
322 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
323
324 assertEquals(0, handler.getLastError().getReal(), epsilonLast);
325 assertEquals(0, handler.getMaximalValueError().getReal(), epsilonMaxValue);
326 assertEquals(0, handler.getMaximalTimeError().getReal(), epsilonMaxTime);
327 assertEquals(name, integ.getName());
328
329 }
330
331 @Test
332 public abstract void testBigStep();
333
334 protected <T extends CalculusFieldElement<T>> void doTestBigStep(Field<T> field,
335 final double belowLast,
336 final double belowMaxValue,
337 final double epsilonMaxTime,
338 final String name)
339 throws MathIllegalArgumentException, MathIllegalStateException {
340
341 TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
342 T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.2);
343
344 FixedStepRungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step);
345 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
346 integ.addStepHandler(handler);
347 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
348
349 assertTrue(handler.getLastError().getReal() > belowLast);
350 assertTrue(handler.getMaximalValueError().getReal() > belowMaxValue);
351 assertEquals(0, handler.getMaximalTimeError().getReal(), epsilonMaxTime);
352 assertEquals(name, integ.getName());
353
354 }
355
356 @Test
357 public abstract void testBackward();
358
359 protected <T extends CalculusFieldElement<T>> void doTestBackward(Field<T> field,
360 final double epsilonLast,
361 final double epsilonMaxValue,
362 final double epsilonMaxTime,
363 final String name)
364 throws MathIllegalArgumentException, MathIllegalStateException {
365
366 TestFieldProblem5<T> pb = new TestFieldProblem5<T>(field);
367 T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.001).abs();
368
369 FixedStepRungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step);
370 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
371 integ.addStepHandler(handler);
372 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
373
374 assertEquals(0, handler.getLastError().getReal(), epsilonLast);
375 assertEquals(0, handler.getMaximalValueError().getReal(), epsilonMaxValue);
376 assertEquals(0, handler.getMaximalTimeError().getReal(), epsilonMaxTime);
377 assertEquals(name, integ.getName());
378
379 }
380
381 @Test
382 public abstract void testKepler();
383
384 protected <T extends CalculusFieldElement<T>> void doTestKepler(Field<T> field, double expectedMaxError, double epsilon)
385 throws MathIllegalArgumentException, MathIllegalStateException {
386
387 final TestFieldProblem3<T> pb = new TestFieldProblem3<T>(field.getZero().add(0.9));
388 T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.0003);
389
390 FixedStepRungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step);
391 integ.addStepHandler(new KeplerHandler<T>(pb, expectedMaxError, epsilon));
392 final FieldExpandableODE<T> expandable = new FieldExpandableODE<T>(pb);
393 assertSame(pb, expandable.getPrimary());
394 integ.integrate(expandable, pb.getInitialState(), pb.getFinalTime());
395 }
396
397 private static class KeplerHandler<T extends CalculusFieldElement<T>> implements FieldODEStepHandler<T> {
398 private T maxError;
399 private final TestFieldProblem3<T> pb;
400 private final double expectedMaxError;
401 private final double epsilon;
402 public KeplerHandler(TestFieldProblem3<T> pb, double expectedMaxError, double epsilon) {
403 this.pb = pb;
404 this.expectedMaxError = expectedMaxError;
405 this.epsilon = epsilon;
406 maxError = pb.getField().getZero();
407 }
408 public void init(FieldODEStateAndDerivative<T> state0, T t) {
409 maxError = pb.getField().getZero();
410 }
411 public void handleStep(FieldODEStateInterpolator<T> interpolator) {
412
413 FieldODEStateAndDerivative<T> current = interpolator.getCurrentState();
414 T[] theoreticalY = pb.computeTheoreticalState(current.getTime());
415 T dx = current.getPrimaryState()[0].subtract(theoreticalY[0]);
416 T dy = current.getPrimaryState()[1].subtract(theoreticalY[1]);
417 T error = dx.square().add(dy.square());
418 if (error.subtract(maxError).getReal() > 0) {
419 maxError = error;
420 }
421 }
422 public void finish(FieldODEStateAndDerivative<T> finalState) {
423 assertEquals(expectedMaxError, maxError.getReal(), epsilon);
424 }
425 }
426
427 @Test
428 public abstract void testStepSize();
429
430 protected <T extends CalculusFieldElement<T>> void doTestStepSize(final Field<T> field, final double epsilon)
431 throws MathIllegalArgumentException, MathIllegalStateException {
432 final T finalTime = field.getZero().add(5.0);
433 final T step = field.getZero().add(1.23456);
434 FixedStepRungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step);
435 integ.addStepHandler(new FieldODEStepHandler<T>() {
436 public void handleStep(FieldODEStateInterpolator<T> interpolator) {
437 if (interpolator.getCurrentState().getTime().subtract(finalTime).getReal() < -0.001) {
438 assertEquals(step.getReal(),
439 interpolator.getCurrentState().getTime().subtract(interpolator.getPreviousState().getTime()).getReal(),
440 epsilon);
441 }
442 }
443 });
444 integ.integrate(new FieldExpandableODE<T>(new FieldOrdinaryDifferentialEquation<T>() {
445 public T[] computeDerivatives(T t, T[] y) {
446 T[] dot = MathArrays.buildArray(t.getField(), 1);
447 dot[0] = t.getField().getOne();
448 return dot;
449 }
450 public int getDimension() {
451 return 1;
452 }
453 }), new FieldODEState<T>(field.getZero(), MathArrays.buildArray(field, 1)), finalTime);
454 }
455
456 @Test
457 public abstract void testSingleStep();
458
459 protected <T extends CalculusFieldElement<T>> void doTestSingleStep(final Field<T> field, final double epsilon) {
460
461 final TestFieldProblem3<T> pb = new TestFieldProblem3<T>(field.getZero().add(0.9));
462 T h = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.0003);
463
464 FixedStepRungeKuttaFieldIntegrator<T> integ = createIntegrator(field, field.getZero().add(Double.NaN));
465 T t = pb.getInitialState().getTime();
466 T[] y = pb.getInitialState().getPrimaryState();
467 for (int i = 0; i < 100; ++i) {
468 y = integ.singleStep(pb, t, y, t.add(h));
469 t = t.add(h);
470 }
471 T[] yth = pb.computeTheoreticalState(t);
472 T dx = y[0].subtract(yth[0]);
473 T dy = y[1].subtract(yth[1]);
474 T error = dx.square().add(dy.square());
475 assertEquals(0.0, error.getReal(), epsilon);
476 }
477
478 @Test
479 public abstract void testTooLargeFirstStep();
480
481 protected <T extends CalculusFieldElement<T>> void doTestTooLargeFirstStep(final Field<T> field) {
482
483 FixedStepRungeKuttaFieldIntegrator<T> integ = createIntegrator(field, field.getZero().add(0.5));
484 final T t0 = field.getZero();
485 final T[] y0 = MathArrays.buildArray(field, 1);
486 y0[0] = field.getOne();
487 final T t = field.getZero().add(0.001);
488 FieldOrdinaryDifferentialEquation<T> equations = new FieldOrdinaryDifferentialEquation<T>() {
489
490 public int getDimension() {
491 return 1;
492 }
493
494 public T[] computeDerivatives(T t, T[] y) {
495 assertTrue(t.getReal() >= FastMath.nextAfter(t0.getReal(), Double.NEGATIVE_INFINITY));
496 assertTrue(t.getReal() <= FastMath.nextAfter(t.getReal(), Double.POSITIVE_INFINITY));
497 T[] yDot = MathArrays.buildArray(field, 1);
498 yDot[0] = y[0].multiply(-100.0);
499 return yDot;
500 }
501
502 };
503
504 integ.integrate(new FieldExpandableODE<T>(equations), new FieldODEState<T>(t0, y0), t);
505
506 }
507
508 @Test
509 public abstract void testUnstableDerivative();
510
511 protected <T extends CalculusFieldElement<T>> void doTestUnstableDerivative(Field<T> field, double epsilon) {
512 final StepFieldProblem<T> stepProblem = new StepFieldProblem<T>(field,
513 (s, isForward) -> 999.0,
514 field.getZero().newInstance(1.0e+12),
515 1000000,
516 field.getZero().newInstance(0.0),
517 field.getZero().newInstance(1.0),
518 field.getZero().newInstance(2.0)).
519 withMaxCheck(field.getZero().newInstance(1.0)).
520 withMaxIter(1000).
521 withThreshold(field.getZero().newInstance(1.0e-12));
522 assertEquals(1.0, stepProblem.getMaxCheckInterval().currentInterval(null, true), 1.0e-15);
523 assertEquals(1000, stepProblem.getMaxIterationCount());
524 assertEquals(1.0e-12, stepProblem.getSolver().getAbsoluteAccuracy().getReal(), 1.0e-25);
525 assertNotNull(stepProblem.getHandler());
526 FixedStepRungeKuttaFieldIntegrator<T> integ = createIntegrator(field, field.getZero().add(0.3));
527 integ.addEventDetector(stepProblem);
528 FieldODEStateAndDerivative<T> result = integ.integrate(new FieldExpandableODE<T>(stepProblem),
529 new FieldODEState<T>(field.getZero(), MathArrays.buildArray(field, 1)),
530 field.getZero().add(10.0));
531 assertEquals(8.0, result.getPrimaryState()[0].getReal(), epsilon);
532 }
533
534 @Test
535 public abstract void testDerivativesConsistency();
536
537 protected <T extends CalculusFieldElement<T>> void doTestDerivativesConsistency(final Field<T> field, double epsilon) {
538 TestFieldProblem3<T> pb = new TestFieldProblem3<T>(field);
539 T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.001);
540 FixedStepRungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step);
541 StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 1.0e-10);
542 }
543
544 @Test
545 public abstract void testPartialDerivatives();
546
547 protected void doTestPartialDerivatives(final double epsilonY,
548 final double[] epsilonPartials) {
549
550
551 final DSFactory factory = new DSFactory(5, 1);
552 final int parOmega = 0;
553 final int parTO = 1;
554 final int parY00 = 2;
555 final int parY01 = 3;
556 final int parT = 4;
557
558 DerivativeStructure omega = factory.variable(parOmega, 1.3);
559 DerivativeStructure t0 = factory.variable(parTO, 1.3);
560 DerivativeStructure[] y0 = new DerivativeStructure[] {
561 factory.variable(parY00, 3.0),
562 factory.variable(parY01, 4.0)
563 };
564 DerivativeStructure t = factory.variable(parT, 6.0);
565 SinCos sinCos = new SinCos(omega);
566
567 FixedStepRungeKuttaFieldIntegrator<DerivativeStructure> integrator =
568 createIntegrator(omega.getField(), t.subtract(t0).multiply(0.001));
569 FieldODEStateAndDerivative<DerivativeStructure> result =
570 integrator.integrate(new FieldExpandableODE<DerivativeStructure>(sinCos),
571 new FieldODEState<DerivativeStructure>(t0, y0),
572 t);
573
574
575 for (int i = 0; i < sinCos.getDimension(); ++i) {
576 assertEquals(sinCos.theoreticalY(t.getReal())[i], result.getPrimaryState()[i].getValue(), epsilonY);
577 }
578
579
580 final double[][] derivatives = sinCos.getDerivatives(t.getReal());
581 for (int i = 0; i < sinCos.getDimension(); ++i) {
582 for (int parameter = 0; parameter < factory.getCompiler().getFreeParameters(); ++parameter) {
583 assertEquals(derivatives[i][parameter],
584 dYdP(result.getPrimaryState()[i], parameter),
585 epsilonPartials[parameter]);
586 }
587 }
588
589 }
590
591 @Test
592 public abstract void testSecondaryEquations();
593
594 protected <T extends CalculusFieldElement<T>> void doTestSecondaryEquations(final Field<T> field,
595 final double epsilonSinCos,
596 final double epsilonLinear) {
597 FieldOrdinaryDifferentialEquation<T> sinCos = new FieldOrdinaryDifferentialEquation<T>() {
598
599 @Override
600 public int getDimension() {
601 return 2;
602 }
603
604 @Override
605 public T[] computeDerivatives(T t, T[] y) {
606
607
608
609
610
611 T[] yDot = y.clone();
612 yDot[0] = y[1].multiply( 0.5);
613 yDot[1] = y[0].multiply(-0.5);
614 return yDot;
615 }
616
617 };
618
619 FieldSecondaryODE<T> linear = new FieldSecondaryODE<T>() {
620
621 @Override
622 public int getDimension() {
623 return 1;
624 }
625
626 @Override
627 public T[] computeDerivatives(T t, T[] primary, T[] primaryDot, T[] secondary) {
628 for (int i = 0; i < primaryDot.length; ++i) {
629
630
631
632
633 primaryDot[i] = primaryDot[i].multiply(2);
634 }
635 T[] secondaryDot = secondary.clone();
636 secondaryDot[0] = t.getField().getOne().negate();
637 return secondaryDot;
638 }
639
640 };
641
642 FieldExpandableODE<T> expandable = new FieldExpandableODE<>(sinCos);
643 expandable.addSecondaryEquations(linear);
644
645 FieldODEIntegrator<T> integrator = createIntegrator(field, field.getZero().add(0.001));
646 final double[] max = new double[2];
647 integrator.addStepHandler(new FieldODEStepHandler<T>() {
648 @Override
649 public void handleStep(FieldODEStateInterpolator<T> interpolator) {
650 for (int i = 0; i <= 10; ++i) {
651 T tPrev = interpolator.getPreviousState().getTime();
652 T tCurr = interpolator.getCurrentState().getTime();
653 T t = tPrev.multiply(10 - i).add(tCurr.multiply(i)).divide(10);
654 FieldODEStateAndDerivative<T> state = interpolator.getInterpolatedState(t);
655 assertEquals(2, state.getPrimaryStateDimension());
656 assertEquals(1, state.getNumberOfSecondaryStates());
657 assertEquals(2, state.getSecondaryStateDimension(0));
658 assertEquals(1, state.getSecondaryStateDimension(1));
659 assertEquals(3, state.getCompleteStateDimension());
660 max[0] = FastMath.max(max[0],
661 t.sin().subtract(state.getPrimaryState()[0]).norm());
662 max[0] = FastMath.max(max[0],
663 t.cos().subtract(state.getPrimaryState()[1]).norm());
664 max[1] = FastMath.max(max[1],
665 field.getOne().subtract(t).subtract(state.getSecondaryState(1)[0]).norm());
666 }
667 }
668 });
669
670 T[] primary0 = MathArrays.buildArray(field, 2);
671 primary0[0] = field.getZero();
672 primary0[1] = field.getOne();
673 T[][] secondary0 = MathArrays.buildArray(field, 1, 1);
674 secondary0[0][0] = field.getOne();
675 FieldODEState<T> initialState = new FieldODEState<T>(field.getZero(), primary0, secondary0);
676
677 FieldODEStateAndDerivative<T> finalState =
678 integrator.integrate(expandable, initialState, field.getZero().add(10.0));
679 assertEquals(10.0, finalState.getTime().getReal(), 1.0e-12);
680 assertEquals(0, max[0], epsilonSinCos);
681 assertEquals(0, max[1], epsilonLinear);
682
683 }
684
685 private double dYdP(final DerivativeStructure y, final int parameter) {
686 int[] orders = new int[y.getFreeParameters()];
687 orders[parameter] = 1;
688 return y.getPartialDerivative(orders);
689 }
690
691 private static class SinCos implements FieldOrdinaryDifferentialEquation<DerivativeStructure> {
692
693 private final DerivativeStructure omega;
694 private DerivativeStructure r;
695 private DerivativeStructure alpha;
696
697 private double dRdY00;
698 private double dRdY01;
699 private double dAlphadOmega;
700 private double dAlphadT0;
701 private double dAlphadY00;
702 private double dAlphadY01;
703
704 protected SinCos(final DerivativeStructure omega) {
705 this.omega = omega;
706 }
707
708 public int getDimension() {
709 return 2;
710 }
711
712 public void init(final DerivativeStructure t0, final DerivativeStructure[] y0,
713 final DerivativeStructure finalTime) {
714
715
716
717 final DerivativeStructure r2 = y0[0].multiply(y0[0]).add(y0[1].multiply(y0[1]));
718
719 this.r = r2.sqrt();
720 this.dRdY00 = y0[0].divide(r).getReal();
721 this.dRdY01 = y0[1].divide(r).getReal();
722
723 this.alpha = y0[0].atan2(y0[1]).subtract(t0.multiply(omega));
724 this.dAlphadOmega = -t0.getReal();
725 this.dAlphadT0 = -omega.getReal();
726 this.dAlphadY00 = y0[1].divide(r2).getReal();
727 this.dAlphadY01 = y0[0].negate().divide(r2).getReal();
728
729 }
730
731 public DerivativeStructure[] computeDerivatives(final DerivativeStructure t, final DerivativeStructure[] y) {
732 return new DerivativeStructure[] {
733 omega.multiply(y[1]),
734 omega.multiply(y[0]).negate()
735 };
736 }
737
738 public double[] theoreticalY(final double t) {
739 final double theta = omega.getReal() * t + alpha.getReal();
740 return new double[] {
741 r.getReal() * FastMath.sin(theta), r.getReal() * FastMath.cos(theta)
742 };
743 }
744
745 public double[][] getDerivatives(final double t) {
746
747
748 final double theta = omega.getReal() * t + alpha.getReal();
749 final double sin = FastMath.sin(theta);
750 final double cos = FastMath.cos(theta);
751 final double y0 = r.getReal() * sin;
752 final double y1 = r.getReal() * cos;
753
754
755 final double dY0dOmega = y1 * (t + dAlphadOmega);
756 final double dY0dT0 = y1 * dAlphadT0;
757 final double dY0dY00 = dRdY00 * sin + y1 * dAlphadY00;
758 final double dY0dY01 = dRdY01 * sin + y1 * dAlphadY01;
759 final double dY0dT = y1 * omega.getReal();
760
761
762 final double dY1dOmega = - y0 * (t + dAlphadOmega);
763 final double dY1dT0 = - y0 * dAlphadT0;
764 final double dY1dY00 = dRdY00 * cos - y0 * dAlphadY00;
765 final double dY1dY01 = dRdY01 * cos - y0 * dAlphadY01;
766 final double dY1dT = - y0 * omega.getReal();
767
768 return new double[][] {
769 { dY0dOmega, dY0dT0, dY0dY00, dY0dY01, dY0dT },
770 { dY1dOmega, dY1dT0, dY1dY00, dY1dY01, dY1dT }
771 };
772
773 }
774
775 }
776
777 @Test
778 public void testIssue250() {
779 final Gradient defaultStep = Gradient.constant(3, 60.);
780 FixedStepRungeKuttaFieldIntegrator<Gradient> integrator = createIntegrator(defaultStep.getField(), defaultStep);
781 assertEquals(defaultStep.getReal(), integrator.getDefaultStep().getReal(), 0.);
782 }
783
784 @Test
785 public void testUsingFieldCoefficients()
786 throws MathIllegalArgumentException, MathIllegalStateException {
787
788 final ComplexField field = ComplexField.getInstance();
789 final TestFieldProblem1<Complex> pb = new TestFieldProblem1<>(field);
790 final double step = FastMath.abs(0.001 * (pb.getFinalTime().getReal() - pb.getInitialState().getTime().getReal()));
791 final Complex fieldStep = new Complex(step);
792
793 final FixedStepRungeKuttaFieldIntegrator<Complex> integratorUsingFieldCoefficients = createIntegrator(field, fieldStep);
794 integratorUsingFieldCoefficients.setUsingFieldCoefficients(true);
795 final FieldODEStateAndDerivative<Complex> terminalState1 = integratorUsingFieldCoefficients.integrate(new FieldExpandableODE<>(pb),
796 pb.getInitialState(), pb.getFinalTime());
797 final FixedStepRungeKuttaFieldIntegrator<Complex> integratorNotUsingFieldCoefficients = createIntegrator(field, fieldStep);
798 integratorNotUsingFieldCoefficients.setUsingFieldCoefficients(false);
799 final FieldODEStateAndDerivative<Complex> terminalState2 = integratorNotUsingFieldCoefficients.integrate(new FieldExpandableODE<>(pb),
800 pb.getInitialState(), pb.getFinalTime());
801
802 final int size = terminalState1.getCompleteStateDimension();
803 for (int i = 0; i < size; i++) {
804 assertEquals(terminalState1.getCompleteState()[i], terminalState2.getCompleteState()[i]);
805 }
806 }
807
808 }