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