1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math4.legacy.ode.nonstiff;
19
20
21 import java.lang.reflect.Array;
22
23 import org.apache.commons.math4.legacy.core.Field;
24 import org.apache.commons.math4.legacy.core.RealFieldElement;
25 import org.apache.commons.math4.legacy.analysis.differentiation.DerivativeStructure;
26 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
27 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
28 import org.apache.commons.math4.legacy.exception.NoBracketingException;
29 import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
30 import org.apache.commons.math4.legacy.ode.FieldExpandableODE;
31 import org.apache.commons.math4.legacy.ode.FirstOrderFieldDifferentialEquations;
32 import org.apache.commons.math4.legacy.ode.FieldODEState;
33 import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative;
34 import org.apache.commons.math4.legacy.ode.TestFieldProblem1;
35 import org.apache.commons.math4.legacy.ode.TestFieldProblem2;
36 import org.apache.commons.math4.legacy.ode.TestFieldProblem3;
37 import org.apache.commons.math4.legacy.ode.TestFieldProblem4;
38 import org.apache.commons.math4.legacy.ode.TestFieldProblem5;
39 import org.apache.commons.math4.legacy.ode.TestFieldProblem6;
40 import org.apache.commons.math4.legacy.ode.TestFieldProblemAbstract;
41 import org.apache.commons.math4.legacy.ode.TestFieldProblemHandler;
42 import org.apache.commons.math4.legacy.ode.events.Action;
43 import org.apache.commons.math4.legacy.ode.events.FieldEventHandler;
44 import org.apache.commons.math4.legacy.ode.sampling.FieldStepHandler;
45 import org.apache.commons.math4.legacy.ode.sampling.FieldStepInterpolator;
46 import org.apache.commons.math4.legacy.ode.sampling.StepInterpolatorTestUtils;
47 import org.apache.commons.math4.core.jdkmath.JdkMath;
48 import org.apache.commons.math4.legacy.core.MathArrays;
49 import org.junit.Assert;
50 import org.junit.Test;
51
52 public abstract class AbstractRungeKuttaFieldIntegratorTest {
53
54 protected abstract <T extends RealFieldElement<T>> RungeKuttaFieldIntegrator<T>
55 createIntegrator(Field<T> field, T step);
56
57 @Test
58 public abstract void testNonFieldIntegratorConsistency();
59
60 protected <T extends RealFieldElement<T>> void doTestNonFieldIntegratorConsistency(final Field<T> field) {
61 try {
62
63
64 RungeKuttaFieldIntegrator<T> fieldIntegrator = createIntegrator(field, field.getZero().add(1));
65 T[][] fieldA = fieldIntegrator.getA();
66 T[] fieldB = fieldIntegrator.getB();
67 T[] fieldC = fieldIntegrator.getC();
68
69 String fieldName = fieldIntegrator.getClass().getName();
70 String regularName = fieldName.replaceAll("Field", "");
71
72
73 @SuppressWarnings("unchecked")
74 Class<RungeKuttaIntegrator> c = (Class<RungeKuttaIntegrator>) Class.forName(regularName);
75 java.lang.reflect.Field jlrFieldA = c.getDeclaredField("STATIC_A");
76 jlrFieldA.setAccessible(true);
77 double[][] regularA = (double[][]) jlrFieldA.get(null);
78 java.lang.reflect.Field jlrFieldB = c.getDeclaredField("STATIC_B");
79 jlrFieldB.setAccessible(true);
80 double[] regularB = (double[]) jlrFieldB.get(null);
81 java.lang.reflect.Field jlrFieldC = c.getDeclaredField("STATIC_C");
82 jlrFieldC.setAccessible(true);
83 double[] regularC = (double[]) jlrFieldC.get(null);
84
85 Assert.assertEquals(regularA.length, fieldA.length);
86 for (int i = 0; i < regularA.length; ++i) {
87 checkArray(regularA[i], fieldA[i]);
88 }
89 checkArray(regularB, fieldB);
90 checkArray(regularC, fieldC);
91 } catch (ClassNotFoundException cnfe) {
92 Assert.fail(cnfe.getLocalizedMessage());
93 } catch (IllegalAccessException iae) {
94 Assert.fail(iae.getLocalizedMessage());
95 } catch (IllegalArgumentException iae) {
96 Assert.fail(iae.getLocalizedMessage());
97 } catch (SecurityException se) {
98 Assert.fail(se.getLocalizedMessage());
99 } catch (NoSuchFieldException nsfe) {
100 Assert.fail(nsfe.getLocalizedMessage());
101 }
102 }
103
104 private <T extends RealFieldElement<T>> void checkArray(double[] regularArray, T[] fieldArray) {
105 Assert.assertEquals(regularArray.length, fieldArray.length);
106 for (int i = 0; i < regularArray.length; ++i) {
107 if (regularArray[i] == 0) {
108 Assert.assertEquals(0.0, fieldArray[i].getReal(), 0.0);
109 } else {
110 Assert.assertEquals(regularArray[i], fieldArray[i].getReal(), JdkMath.ulp(regularArray[i]));
111 }
112 }
113 }
114
115 @Test
116 public abstract void testMissedEndEvent();
117
118 protected <T extends RealFieldElement<T>> void doTestMissedEndEvent(final Field<T> field,
119 final double epsilonT, final double epsilonY)
120 throws DimensionMismatchException, NumberIsTooSmallException,
121 MaxCountExceededException, NoBracketingException {
122 final T t0 = field.getZero().add(1878250320.0000029);
123 final T tEvent = field.getZero().add(1878250379.9999986);
124 final T[] k = MathArrays.buildArray(field, 3);
125 k[0] = field.getZero().add(1.0e-4);
126 k[1] = field.getZero().add(1.0e-5);
127 k[2] = field.getZero().add(1.0e-6);
128 FirstOrderFieldDifferentialEquations<T> ode = new FirstOrderFieldDifferentialEquations<T>() {
129
130 @Override
131 public int getDimension() {
132 return k.length;
133 }
134
135 @Override
136 public void init(T t0, T[] y0, T t) {
137 }
138
139 @Override
140 public T[] computeDerivatives(T t, T[] y) {
141 T[] yDot = MathArrays.buildArray(field, k.length);
142 for (int i = 0; i < y.length; ++i) {
143 yDot[i] = k[i].multiply(y[i]);
144 }
145 return yDot;
146 }
147 };
148
149 RungeKuttaFieldIntegrator<T> integrator = createIntegrator(field, field.getZero().add(60.0));
150
151 T[] y0 = MathArrays.buildArray(field, k.length);
152 for (int i = 0; i < y0.length; ++i) {
153 y0[i] = field.getOne().add(i);
154 }
155
156 FieldODEStateAndDerivative<T> result = integrator.integrate(new FieldExpandableODE<>(ode),
157 new FieldODEState<>(t0, y0),
158 tEvent);
159 Assert.assertEquals(tEvent.getReal(), result.getTime().getReal(), epsilonT);
160 T[] y = result.getState();
161 for (int i = 0; i < y.length; ++i) {
162 Assert.assertEquals(y0[i].multiply(k[i].multiply(result.getTime().subtract(t0)).exp()).getReal(),
163 y[i].getReal(),
164 epsilonY);
165 }
166
167 integrator.addEventHandler(new FieldEventHandler<T>() {
168
169 @Override
170 public void init(FieldODEStateAndDerivative<T> state0, T t) {
171 }
172
173 @Override
174 public FieldODEState<T> resetState(FieldODEStateAndDerivative<T> state) {
175 return state;
176 }
177
178 @Override
179 public T g(FieldODEStateAndDerivative<T> state) {
180 return state.getTime().subtract(tEvent);
181 }
182
183 @Override
184 public Action eventOccurred(FieldODEStateAndDerivative<T> state, boolean increasing) {
185 Assert.assertEquals(tEvent.getReal(), state.getTime().getReal(), epsilonT);
186 return Action.CONTINUE;
187 }
188 }, Double.POSITIVE_INFINITY, 1.0e-20, 100);
189 result = integrator.integrate(new FieldExpandableODE<>(ode),
190 new FieldODEState<>(t0, y0),
191 tEvent.add(120));
192 Assert.assertEquals(tEvent.add(120).getReal(), result.getTime().getReal(), epsilonT);
193 y = result.getState();
194 for (int i = 0; i < y.length; ++i) {
195 Assert.assertEquals(y0[i].multiply(k[i].multiply(result.getTime().subtract(t0)).exp()).getReal(),
196 y[i].getReal(),
197 epsilonY);
198 }
199 }
200
201 @Test
202 public abstract void testSanityChecks();
203
204 protected <T extends RealFieldElement<T>> void doTestSanityChecks(Field<T> field)
205 throws DimensionMismatchException, NumberIsTooSmallException,
206 MaxCountExceededException, NoBracketingException {
207 RungeKuttaFieldIntegrator<T> integrator = createIntegrator(field, field.getZero().add(0.01));
208 try {
209 TestFieldProblem1<T> pb = new TestFieldProblem1<>(field);
210 integrator.integrate(new FieldExpandableODE<>(pb),
211 new FieldODEState<>(field.getZero(), MathArrays.buildArray(field, pb.getDimension() + 10)),
212 field.getOne());
213 Assert.fail("an exception should have been thrown");
214 } catch(DimensionMismatchException ie) {
215 }
216 try {
217 TestFieldProblem1<T> pb = new TestFieldProblem1<>(field);
218 integrator.integrate(new FieldExpandableODE<>(pb),
219 new FieldODEState<>(field.getZero(), MathArrays.buildArray(field, pb.getDimension())),
220 field.getZero());
221 Assert.fail("an exception should have been thrown");
222 } catch(NumberIsTooSmallException ie) {
223 }
224 }
225
226 @Test
227 public abstract void testDecreasingSteps();
228
229 protected <T extends RealFieldElement<T>> void doTestDecreasingSteps(Field<T> field,
230 final double safetyValueFactor,
231 final double safetyTimeFactor,
232 final double epsilonT)
233 throws DimensionMismatchException, NumberIsTooSmallException,
234 MaxCountExceededException, NoBracketingException {
235
236 @SuppressWarnings("unchecked")
237 TestFieldProblemAbstract<T>[] allProblems =
238 (TestFieldProblemAbstract<T>[]) Array.newInstance(TestFieldProblemAbstract.class, 6);
239 allProblems[0] = new TestFieldProblem1<>(field);
240 allProblems[1] = new TestFieldProblem2<>(field);
241 allProblems[2] = new TestFieldProblem3<>(field);
242 allProblems[3] = new TestFieldProblem4<>(field);
243 allProblems[4] = new TestFieldProblem5<>(field);
244 allProblems[5] = new TestFieldProblem6<>(field);
245 for (TestFieldProblemAbstract<T> pb : allProblems) {
246
247 T previousValueError = null;
248 T previousTimeError = null;
249 for (int i = 4; i < 10; ++i) {
250
251 T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(JdkMath.pow(2.0, -i));
252
253 RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step);
254 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<>(pb, integ);
255 integ.addStepHandler(handler);
256 FieldEventHandler<T>[] functions = pb.getEventsHandlers();
257 for (int l = 0; l < functions.length; ++l) {
258 integ.addEventHandler(functions[l],
259 Double.POSITIVE_INFINITY, 1.0e-6 * step.getReal(), 1000);
260 }
261 Assert.assertEquals(functions.length, integ.getEventHandlers().size());
262 FieldODEStateAndDerivative<T> stop = integ.integrate(new FieldExpandableODE<>(pb),
263 pb.getInitialState(),
264 pb.getFinalTime());
265 if (functions.length == 0) {
266 Assert.assertEquals(pb.getFinalTime().getReal(), stop.getTime().getReal(), epsilonT);
267 }
268
269 T error = handler.getMaximalValueError();
270 if (i > 4) {
271 Assert.assertTrue(error.subtract(previousValueError.abs().multiply(safetyValueFactor)).getReal() < 0);
272 }
273 previousValueError = error;
274
275 T timeError = handler.getMaximalTimeError();
276 if (i > 4) {
277 Assert.assertTrue(timeError.subtract(previousTimeError.abs().multiply(safetyTimeFactor)).getReal() <= 0);
278 }
279 previousTimeError = timeError;
280
281 integ.clearEventHandlers();
282 Assert.assertEquals(0, integ.getEventHandlers().size());
283 }
284 }
285 }
286
287 @Test
288 public abstract void testSmallStep();
289
290 protected <T extends RealFieldElement<T>> void doTestSmallStep(Field<T> field,
291 final double epsilonLast,
292 final double epsilonMaxValue,
293 final double epsilonMaxTime,
294 final String name)
295 throws DimensionMismatchException, NumberIsTooSmallException,
296 MaxCountExceededException, NoBracketingException {
297
298 TestFieldProblem1<T> pb = new TestFieldProblem1<>(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<>(pb, integ);
303 integ.addStepHandler(handler);
304 integ.integrate(new FieldExpandableODE<>(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 @Test
313 public abstract void testBigStep();
314
315 protected <T extends RealFieldElement<T>> void doTestBigStep(Field<T> field,
316 final double belowLast,
317 final double belowMaxValue,
318 final double epsilonMaxTime,
319 final String name)
320 throws DimensionMismatchException, NumberIsTooSmallException,
321 MaxCountExceededException, NoBracketingException {
322
323 TestFieldProblem1<T> pb = new TestFieldProblem1<>(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<>(pb, integ);
328 integ.addStepHandler(handler);
329 integ.integrate(new FieldExpandableODE<>(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 @Test
338 public abstract void testBackward();
339
340 protected <T extends RealFieldElement<T>> void doTestBackward(Field<T> field,
341 final double epsilonLast,
342 final double epsilonMaxValue,
343 final double epsilonMaxTime,
344 final String name)
345 throws DimensionMismatchException, NumberIsTooSmallException,
346 MaxCountExceededException, NoBracketingException {
347
348 TestFieldProblem5<T> pb = new TestFieldProblem5<>(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<>(pb, integ);
353 integ.addStepHandler(handler);
354 integ.integrate(new FieldExpandableODE<>(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 @Test
363 public abstract void testKepler();
364
365 protected <T extends RealFieldElement<T>> void doTestKepler(Field<T> field, double expectedMaxError, double epsilon)
366 throws DimensionMismatchException, NumberIsTooSmallException,
367 MaxCountExceededException, NoBracketingException {
368
369 final TestFieldProblem3<T> pb = new TestFieldProblem3<>(field, 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<>(pb, expectedMaxError, epsilon));
374 integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
375 }
376
377 private static final class KeplerHandler<T extends RealFieldElement<T>> implements FieldStepHandler<T> {
378 private T maxError;
379 private final TestFieldProblem3<T> pb;
380 private final double expectedMaxError;
381 private final double epsilon;
382 KeplerHandler(TestFieldProblem3<T> pb, double expectedMaxError, double epsilon) {
383 this.pb = pb;
384 this.expectedMaxError = expectedMaxError;
385 this.epsilon = epsilon;
386 maxError = pb.getField().getZero();
387 }
388 @Override
389 public void init(FieldODEStateAndDerivative<T> state0, T t) {
390 maxError = pb.getField().getZero();
391 }
392 @Override
393 public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast)
394 throws MaxCountExceededException {
395
396 FieldODEStateAndDerivative<T> current = interpolator.getCurrentState();
397 T[] theoreticalY = pb.computeTheoreticalState(current.getTime());
398 T dx = current.getState()[0].subtract(theoreticalY[0]);
399 T dy = current.getState()[1].subtract(theoreticalY[1]);
400 T error = dx.multiply(dx).add(dy.multiply(dy));
401 if (error.subtract(maxError).getReal() > 0) {
402 maxError = error;
403 }
404 if (isLast) {
405 Assert.assertEquals(expectedMaxError, maxError.getReal(), epsilon);
406 }
407 }
408 }
409
410 @Test
411 public abstract void testStepSize();
412
413 protected <T extends RealFieldElement<T>> void doTestStepSize(final Field<T> field, final double epsilon)
414 throws DimensionMismatchException, NumberIsTooSmallException,
415 MaxCountExceededException, NoBracketingException {
416 final T step = field.getZero().add(1.23456);
417 RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step);
418 integ.addStepHandler(new FieldStepHandler<T>() {
419 @Override
420 public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast) {
421 if (! isLast) {
422 Assert.assertEquals(step.getReal(),
423 interpolator.getCurrentState().getTime().subtract(interpolator.getPreviousState().getTime()).getReal(),
424 epsilon);
425 }
426 }
427 @Override
428 public void init(FieldODEStateAndDerivative<T> s0, T t) {
429 }
430 });
431 integ.integrate(new FieldExpandableODE<>(new FirstOrderFieldDifferentialEquations<T>() {
432 @Override
433 public void init(T t0, T[] y0, T t) {
434 }
435 @Override
436 public T[] computeDerivatives(T t, T[] y) {
437 T[] dot = MathArrays.buildArray(t.getField(), 1);
438 dot[0] = t.getField().getOne();
439 return dot;
440 }
441 @Override
442 public int getDimension() {
443 return 1;
444 }
445 }), new FieldODEState<>(field.getZero(), MathArrays.buildArray(field, 1)), field.getZero().add(5.0));
446 }
447
448 @Test
449 public abstract void testSingleStep();
450
451 protected <T extends RealFieldElement<T>> void doTestSingleStep(final Field<T> field, final double epsilon) {
452
453 final TestFieldProblem3<T> pb = new TestFieldProblem3<>(field, field.getZero().add(0.9));
454 T h = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.0003);
455
456 RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, field.getZero().add(Double.NaN));
457 T t = pb.getInitialState().getTime();
458 T[] y = pb.getInitialState().getState();
459 for (int i = 0; i < 100; ++i) {
460 y = integ.singleStep(pb, t, y, t.add(h));
461 t = t.add(h);
462 }
463 T[] yth = pb.computeTheoreticalState(t);
464 T dx = y[0].subtract(yth[0]);
465 T dy = y[1].subtract(yth[1]);
466 T error = dx.multiply(dx).add(dy.multiply(dy));
467 Assert.assertEquals(0.0, error.getReal(), epsilon);
468 }
469
470 @Test
471 public abstract void testTooLargeFirstStep();
472
473 protected <T extends RealFieldElement<T>> void doTestTooLargeFirstStep(final Field<T> field) {
474
475 RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, field.getZero().add(0.5));
476 final T t0 = field.getZero();
477 final T[] y0 = MathArrays.buildArray(field, 1);
478 y0[0] = field.getOne();
479 final T t = field.getZero().add(0.001);
480 FirstOrderFieldDifferentialEquations<T> equations = new FirstOrderFieldDifferentialEquations<T>() {
481
482 @Override
483 public int getDimension() {
484 return 1;
485 }
486
487 @Override
488 public void init(T t0, T[] y0, T t) {
489 }
490
491 @Override
492 public T[] computeDerivatives(T t, T[] y) {
493 Assert.assertTrue(t.getReal() >= JdkMath.nextAfter(t0.getReal(), Double.NEGATIVE_INFINITY));
494 Assert.assertTrue(t.getReal() <= JdkMath.nextAfter(t.getReal(), Double.POSITIVE_INFINITY));
495 T[] yDot = MathArrays.buildArray(field, 1);
496 yDot[0] = y[0].multiply(-100.0);
497 return yDot;
498 }
499 };
500
501 integ.integrate(new FieldExpandableODE<>(equations), new FieldODEState<>(t0, y0), t);
502 }
503
504 @Test
505 public abstract void testUnstableDerivative();
506
507 protected <T extends RealFieldElement<T>> void doTestUnstableDerivative(Field<T> field, double epsilon) {
508 final StepFieldProblem<T> stepProblem = new StepFieldProblem<>(field,
509 field.getZero().add(0.0),
510 field.getZero().add(1.0),
511 field.getZero().add(2.0));
512 RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, field.getZero().add(0.3));
513 integ.addEventHandler(stepProblem, 1.0, 1.0e-12, 1000);
514 FieldODEStateAndDerivative<T> result = integ.integrate(new FieldExpandableODE<>(stepProblem),
515 new FieldODEState<>(field.getZero(), MathArrays.buildArray(field, 1)),
516 field.getZero().add(10.0));
517 Assert.assertEquals(8.0, result.getState()[0].getReal(), epsilon);
518 }
519
520 @Test
521 public abstract void testDerivativesConsistency();
522
523 protected <T extends RealFieldElement<T>> void doTestDerivativesConsistency(final Field<T> field, double epsilon) {
524 TestFieldProblem3<T> pb = new TestFieldProblem3<>(field);
525 T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.001);
526 RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step);
527 StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 1.0e-10);
528 }
529
530 @Test
531 public abstract void testPartialDerivatives();
532
533 protected <T extends RealFieldElement<T>> void doTestPartialDerivatives(final double epsilonY,
534 final double[] epsilonPartials) {
535
536
537 final int parameters = 5;
538 final int order = 1;
539 final int parOmega = 0;
540 final int parTO = 1;
541 final int parY00 = 2;
542 final int parY01 = 3;
543 final int parT = 4;
544
545 DerivativeStructure omega = new DerivativeStructure(parameters, order, parOmega, 1.3);
546 DerivativeStructure t0 = new DerivativeStructure(parameters, order, parTO, 1.3);
547 DerivativeStructure[] y0 = new DerivativeStructure[] {
548 new DerivativeStructure(parameters, order, parY00, 3.0),
549 new DerivativeStructure(parameters, order, parY01, 4.0)
550 };
551 DerivativeStructure t = new DerivativeStructure(parameters, order, parT, 6.0);
552 SinCos sinCos = new SinCos(omega);
553
554 RungeKuttaFieldIntegrator<DerivativeStructure> integrator =
555 createIntegrator(omega.getField(), t.subtract(t0).multiply(0.001));
556 FieldODEStateAndDerivative<DerivativeStructure> result =
557 integrator.integrate(new FieldExpandableODE<>(sinCos),
558 new FieldODEState<>(t0, y0),
559 t);
560
561
562 for (int i = 0; i < sinCos.getDimension(); ++i) {
563 Assert.assertEquals(sinCos.theoreticalY(t.getReal())[i], result.getState()[i].getValue(), epsilonY);
564 }
565
566
567 final double[][] derivatives = sinCos.getDerivatives(t.getReal());
568 for (int i = 0; i < sinCos.getDimension(); ++i) {
569 for (int parameter = 0; parameter < parameters; ++parameter) {
570 Assert.assertEquals(derivatives[i][parameter],
571 dYdP(result.getState()[i], parameter),
572 epsilonPartials[parameter]);
573 }
574 }
575 }
576
577 private double dYdP(final DerivativeStructure y, final int parameter) {
578 int[] orders = new int[y.getFreeParameters()];
579 orders[parameter] = 1;
580 return y.getPartialDerivative(orders);
581 }
582
583 private static final class SinCos implements FirstOrderFieldDifferentialEquations<DerivativeStructure> {
584
585 private final DerivativeStructure omega;
586 private DerivativeStructure r;
587 private DerivativeStructure alpha;
588
589 private double dRdY00;
590 private double dRdY01;
591 private double dAlphadOmega;
592 private double dAlphadT0;
593 private double dAlphadY00;
594 private double dAlphadY01;
595
596 protected SinCos(final DerivativeStructure omega) {
597 this.omega = omega;
598 }
599
600 @Override
601 public int getDimension() {
602 return 2;
603 }
604
605 @Override
606 public void init(final DerivativeStructure t0, final DerivativeStructure[] y0,
607 final DerivativeStructure finalTime) {
608
609
610
611 final DerivativeStructure r2 = y0[0].multiply(y0[0]).add(y0[1].multiply(y0[1]));
612
613 this.r = r2.sqrt();
614 this.dRdY00 = y0[0].divide(r).getReal();
615 this.dRdY01 = y0[1].divide(r).getReal();
616
617 this.alpha = y0[0].atan2(y0[1]).subtract(t0.multiply(omega));
618 this.dAlphadOmega = -t0.getReal();
619 this.dAlphadT0 = -omega.getReal();
620 this.dAlphadY00 = y0[1].divide(r2).getReal();
621 this.dAlphadY01 = y0[0].negate().divide(r2).getReal();
622 }
623
624 @Override
625 public DerivativeStructure[] computeDerivatives(final DerivativeStructure t, final DerivativeStructure[] y) {
626 return new DerivativeStructure[] {
627 omega.multiply(y[1]),
628 omega.multiply(y[0]).negate()
629 };
630 }
631
632 public double[] theoreticalY(final double t) {
633 final double theta = omega.getReal() * t + alpha.getReal();
634 return new double[] {
635 r.getReal() * JdkMath.sin(theta), r.getReal() * JdkMath.cos(theta)
636 };
637 }
638
639 public double[][] getDerivatives(final double t) {
640
641
642 final double theta = omega.getReal() * t + alpha.getReal();
643 final double sin = JdkMath.sin(theta);
644 final double cos = JdkMath.cos(theta);
645 final double y0 = r.getReal() * sin;
646 final double y1 = r.getReal() * cos;
647
648
649 final double dY0dOmega = y1 * (t + dAlphadOmega);
650 final double dY0dT0 = y1 * dAlphadT0;
651 final double dY0dY00 = dRdY00 * sin + y1 * dAlphadY00;
652 final double dY0dY01 = dRdY01 * sin + y1 * dAlphadY01;
653 final double dY0dT = y1 * omega.getReal();
654
655
656 final double dY1dOmega = - y0 * (t + dAlphadOmega);
657 final double dY1dT0 = - y0 * dAlphadT0;
658 final double dY1dY00 = dRdY00 * cos - y0 * dAlphadY00;
659 final double dY1dY01 = dRdY01 * cos - y0 * dAlphadY01;
660 final double dY1dT = - y0 * omega.getReal();
661
662 return new double[][] {
663 { dY0dOmega, dY0dT0, dY0dY00, dY0dY01, dY0dT },
664 { dY1dOmega, dY1dT0, dY1dY00, dY1dY01, dY1dT }
665 };
666 }
667 }
668 }