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 org.apache.commons.math4.legacy.core.Field;
22 import org.apache.commons.math4.legacy.core.RealFieldElement;
23 import org.apache.commons.math4.legacy.analysis.differentiation.DerivativeStructure;
24 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
25 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
26 import org.apache.commons.math4.legacy.exception.NoBracketingException;
27 import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
28 import org.apache.commons.math4.legacy.ode.FieldExpandableODE;
29 import org.apache.commons.math4.legacy.ode.FirstOrderFieldDifferentialEquations;
30 import org.apache.commons.math4.legacy.ode.FirstOrderFieldIntegrator;
31 import org.apache.commons.math4.legacy.ode.FieldODEState;
32 import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative;
33 import org.apache.commons.math4.legacy.ode.TestFieldProblem1;
34 import org.apache.commons.math4.legacy.ode.TestFieldProblem3;
35 import org.apache.commons.math4.legacy.ode.TestFieldProblem4;
36 import org.apache.commons.math4.legacy.ode.TestFieldProblem5;
37 import org.apache.commons.math4.legacy.ode.TestFieldProblemHandler;
38 import org.apache.commons.math4.legacy.ode.events.Action;
39 import org.apache.commons.math4.legacy.ode.events.FieldEventHandler;
40 import org.apache.commons.math4.legacy.ode.sampling.FieldStepHandler;
41 import org.apache.commons.math4.legacy.ode.sampling.FieldStepInterpolator;
42 import org.apache.commons.math4.core.jdkmath.JdkMath;
43 import org.apache.commons.math4.legacy.core.MathArrays;
44 import org.junit.Assert;
45 import org.junit.Test;
46
47 public abstract class AbstractEmbeddedRungeKuttaFieldIntegratorTest {
48
49 protected abstract <T extends RealFieldElement<T>> EmbeddedRungeKuttaFieldIntegrator<T>
50 createIntegrator(Field<T> field, double minStep, double maxStep,
51 double scalAbsoluteTolerance, double scalRelativeTolerance);
52
53 protected abstract <T extends RealFieldElement<T>> EmbeddedRungeKuttaFieldIntegrator<T>
54 createIntegrator(Field<T> field, double minStep, double maxStep,
55 double[] vecAbsoluteTolerance, double[] vecRelativeTolerance);
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 EmbeddedRungeKuttaFieldIntegrator<T> fieldIntegrator = createIntegrator(field, 0.001, 1.0, 1.0, 1.0);
65 T[][] fieldA = fieldIntegrator.getA();
66 T[] fieldB = fieldIntegrator.getB();
67 T[] fieldC = fieldIntegrator.getC();
68 if (fieldIntegrator instanceof DormandPrince853FieldIntegrator) {
69
70
71
72 T[][] reducedFieldA = MathArrays.buildArray(field, 12, -1);
73 T[] reducedFieldB = MathArrays.buildArray(field, 13);
74 T[] reducedFieldC = MathArrays.buildArray(field, 12);
75 System.arraycopy(fieldA, 0, reducedFieldA, 0, reducedFieldA.length);
76 System.arraycopy(fieldB, 0, reducedFieldB, 0, reducedFieldB.length);
77 System.arraycopy(fieldC, 0, reducedFieldC, 0, reducedFieldC.length);
78 fieldA = reducedFieldA;
79 fieldB = reducedFieldB;
80 fieldC = reducedFieldC;
81 }
82
83 String fieldName = fieldIntegrator.getClass().getName();
84 String regularName = fieldName.replaceAll("Field", "");
85
86
87 @SuppressWarnings("unchecked")
88 Class<RungeKuttaIntegrator> c = (Class<RungeKuttaIntegrator>) Class.forName(regularName);
89 java.lang.reflect.Field jlrFieldA = c.getDeclaredField("STATIC_A");
90 jlrFieldA.setAccessible(true);
91 double[][] regularA = (double[][]) jlrFieldA.get(null);
92 java.lang.reflect.Field jlrFieldB = c.getDeclaredField("STATIC_B");
93 jlrFieldB.setAccessible(true);
94 double[] regularB = (double[]) jlrFieldB.get(null);
95 java.lang.reflect.Field jlrFieldC = c.getDeclaredField("STATIC_C");
96 jlrFieldC.setAccessible(true);
97 double[] regularC = (double[]) jlrFieldC.get(null);
98
99 Assert.assertEquals(regularA.length, fieldA.length);
100 for (int i = 0; i < regularA.length; ++i) {
101 checkArray(regularA[i], fieldA[i]);
102 }
103 checkArray(regularB, fieldB);
104 checkArray(regularC, fieldC);
105 } catch (ClassNotFoundException cnfe) {
106 Assert.fail(cnfe.getLocalizedMessage());
107 } catch (IllegalAccessException iae) {
108 Assert.fail(iae.getLocalizedMessage());
109 } catch (IllegalArgumentException iae) {
110 Assert.fail(iae.getLocalizedMessage());
111 } catch (SecurityException se) {
112 Assert.fail(se.getLocalizedMessage());
113 } catch (NoSuchFieldException nsfe) {
114 Assert.fail(nsfe.getLocalizedMessage());
115 }
116 }
117
118 private <T extends RealFieldElement<T>> void checkArray(double[] regularArray, T[] fieldArray) {
119 Assert.assertEquals(regularArray.length, fieldArray.length);
120 for (int i = 0; i < regularArray.length; ++i) {
121 if (regularArray[i] == 0) {
122 Assert.assertEquals(0.0, fieldArray[i].getReal(), 0.0);
123 } else {
124 Assert.assertEquals(regularArray[i], fieldArray[i].getReal(), JdkMath.ulp(regularArray[i]));
125 }
126 }
127 }
128
129 @Test
130 public abstract void testForwardBackwardExceptions();
131
132 protected <T extends RealFieldElement<T>> void doTestForwardBackwardExceptions(final Field<T> field) {
133 FirstOrderFieldDifferentialEquations<T> equations = new FirstOrderFieldDifferentialEquations<T>() {
134
135 @Override
136 public int getDimension() {
137 return 1;
138 }
139
140 @Override
141 public void init(T t0, T[] y0, T t) {
142 }
143
144 @Override
145 public T[] computeDerivatives(T t, T[] y) {
146 if (t.getReal() < -0.5) {
147 throw new LocalException();
148 } else {
149 throw new RuntimeException("oops");
150 }
151 }
152 };
153
154 EmbeddedRungeKuttaFieldIntegrator<T> integrator = createIntegrator(field, 0.0, 1.0, 1.0e-10, 1.0e-10);
155
156 try {
157 integrator.integrate(new FieldExpandableODE<>(equations),
158 new FieldODEState<>(field.getOne().negate(),
159 MathArrays.buildArray(field, 1)),
160 field.getZero());
161 Assert.fail("an exception should have been thrown");
162 } catch(LocalException de) {
163
164 }
165
166 try {
167 integrator.integrate(new FieldExpandableODE<>(equations),
168 new FieldODEState<>(field.getZero(),
169 MathArrays.buildArray(field, 1)),
170 field.getOne());
171 Assert.fail("an exception should have been thrown");
172 } catch(RuntimeException de) {
173
174 }
175 }
176
177 protected static class LocalException extends RuntimeException {
178 private static final long serialVersionUID = 20151208L;
179 }
180
181 @Test(expected=NumberIsTooSmallException.class)
182 public abstract void testMinStep();
183
184 protected <T extends RealFieldElement<T>> void doTestMinStep(final Field<T> field)
185 throws NumberIsTooSmallException {
186
187 TestFieldProblem1<T> pb = new TestFieldProblem1<>(field);
188 double minStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.1).getReal();
189 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
190 double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
191 double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
192
193 FirstOrderFieldIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
194 vecAbsoluteTolerance, vecRelativeTolerance);
195 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<>(pb, integ);
196 integ.addStepHandler(handler);
197 integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
198 Assert.fail("an exception should have been thrown");
199 }
200
201 @Test
202 public abstract void testIncreasingTolerance();
203
204 protected <T extends RealFieldElement<T>> void doTestIncreasingTolerance(final Field<T> field,
205 double factor,
206 double epsilon) {
207
208 int previousCalls = Integer.MAX_VALUE;
209 for (int i = -12; i < -2; ++i) {
210 TestFieldProblem1<T> pb = new TestFieldProblem1<>(field);
211 double minStep = 0;
212 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
213 double scalAbsoluteTolerance = JdkMath.pow(10.0, i);
214 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
215
216 FirstOrderFieldIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
217 scalAbsoluteTolerance, scalRelativeTolerance);
218 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<>(pb, integ);
219 integ.addStepHandler(handler);
220 integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
221
222 Assert.assertTrue(handler.getMaximalValueError().getReal() < (factor * scalAbsoluteTolerance));
223 Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), epsilon);
224
225 int calls = pb.getCalls();
226 Assert.assertEquals(integ.getEvaluations(), calls);
227 Assert.assertTrue(calls <= previousCalls);
228 previousCalls = calls;
229 }
230 }
231
232 @Test
233 public abstract void testEvents();
234
235 protected <T extends RealFieldElement<T>> void doTestEvents(final Field<T> field,
236 final double epsilonMaxValue,
237 final String name) {
238
239 TestFieldProblem4<T> pb = new TestFieldProblem4<>(field);
240 double minStep = 0;
241 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
242 double scalAbsoluteTolerance = 1.0e-8;
243 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
244
245 FirstOrderFieldIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
246 scalAbsoluteTolerance, scalRelativeTolerance);
247 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<>(pb, integ);
248 integ.addStepHandler(handler);
249 FieldEventHandler<T>[] functions = pb.getEventsHandlers();
250 double convergence = 1.0e-8 * maxStep;
251 for (int l = 0; l < functions.length; ++l) {
252 integ.addEventHandler(functions[l], Double.POSITIVE_INFINITY, convergence, 1000);
253 }
254 Assert.assertEquals(functions.length, integ.getEventHandlers().size());
255 integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
256
257 Assert.assertEquals(0, handler.getMaximalValueError().getReal(), epsilonMaxValue);
258 Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), convergence);
259 Assert.assertEquals(12.0, handler.getLastTime().getReal(), convergence);
260 Assert.assertEquals(name, integ.getName());
261 integ.clearEventHandlers();
262 Assert.assertEquals(0, integ.getEventHandlers().size());
263 }
264
265 @Test(expected=LocalException.class)
266 public abstract void testEventsErrors();
267
268 protected <T extends RealFieldElement<T>> void doTestEventsErrors(final Field<T> field)
269 throws LocalException {
270 final TestFieldProblem1<T> pb = new TestFieldProblem1<>(field);
271 double minStep = 0;
272 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
273 double scalAbsoluteTolerance = 1.0e-8;
274 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
275
276 FirstOrderFieldIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
277 scalAbsoluteTolerance, scalRelativeTolerance);
278 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<>(pb, integ);
279 integ.addStepHandler(handler);
280
281 integ.addEventHandler(new FieldEventHandler<T>() {
282 @Override
283 public void init(FieldODEStateAndDerivative<T> state0, T t) {
284 }
285 @Override
286 public Action eventOccurred(FieldODEStateAndDerivative<T> state, boolean increasing) {
287 return Action.CONTINUE;
288 }
289 @Override
290 public T g(FieldODEStateAndDerivative<T> state) {
291 T middle = pb.getInitialState().getTime().add(pb.getFinalTime()).multiply(0.5);
292 T offset = state.getTime().subtract(middle);
293 if (offset.getReal() > 0) {
294 throw new LocalException();
295 }
296 return offset;
297 }
298 @Override
299 public FieldODEState<T> resetState(FieldODEStateAndDerivative<T> state) {
300 return state;
301 }
302 }, Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 1000);
303
304 integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
305 }
306
307 @Test
308 public abstract void testEventsNoConvergence();
309
310 protected <T extends RealFieldElement<T>> void doTestEventsNoConvergence(final Field<T> field){
311
312 final TestFieldProblem1<T> pb = new TestFieldProblem1<>(field);
313 double minStep = 0;
314 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
315 double scalAbsoluteTolerance = 1.0e-8;
316 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
317
318 FirstOrderFieldIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
319 scalAbsoluteTolerance, scalRelativeTolerance);
320 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<>(pb, integ);
321 integ.addStepHandler(handler);
322
323 integ.addEventHandler(new FieldEventHandler<T>() {
324 @Override
325 public void init(FieldODEStateAndDerivative<T> state0, T t) {
326 }
327 @Override
328 public Action eventOccurred(FieldODEStateAndDerivative<T> state, boolean increasing) {
329 return Action.CONTINUE;
330 }
331 @Override
332 public T g(FieldODEStateAndDerivative<T> state) {
333 T middle = pb.getInitialState().getTime().add(pb.getFinalTime()).multiply(0.5);
334 T offset = state.getTime().subtract(middle);
335 return (offset.getReal() > 0) ? offset.add(0.5) : offset.subtract(0.5);
336 }
337 @Override
338 public FieldODEState<T> resetState(FieldODEStateAndDerivative<T> state) {
339 return state;
340 }
341 }, Double.POSITIVE_INFINITY, 1.0e-8 * maxStep, 3);
342
343 try {
344 integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
345 Assert.fail("an exception should have been thrown");
346 } catch (MaxCountExceededException mcee) {
347
348 }
349 }
350
351 @Test
352 public abstract void testSanityChecks();
353
354 protected <T extends RealFieldElement<T>> void doTestSanityChecks(Field<T> field) {
355 TestFieldProblem3<T> pb = new TestFieldProblem3<>(field);
356 try {
357 EmbeddedRungeKuttaFieldIntegrator<T> integrator = createIntegrator(field, 0,
358 pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal(),
359 new double[4], new double[4]);
360 integrator.integrate(new FieldExpandableODE<>(pb),
361 new FieldODEState<>(pb.getInitialState().getTime(),
362 MathArrays.buildArray(field, 6)),
363 pb.getFinalTime());
364 Assert.fail("an exception should have been thrown");
365 } catch(DimensionMismatchException ie) {
366 }
367 try {
368 EmbeddedRungeKuttaFieldIntegrator<T> integrator =
369 createIntegrator(field, 0,
370 pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal(),
371 new double[2], new double[4]);
372 integrator.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
373 Assert.fail("an exception should have been thrown");
374 } catch(DimensionMismatchException ie) {
375 }
376 try {
377 EmbeddedRungeKuttaFieldIntegrator<T> integrator =
378 createIntegrator(field, 0,
379 pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal(),
380 new double[4], new double[4]);
381 integrator.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getInitialState().getTime());
382 Assert.fail("an exception should have been thrown");
383 } catch(NumberIsTooSmallException ie) {
384 }
385 }
386
387 @Test
388 public abstract void testBackward();
389
390 protected <T extends RealFieldElement<T>> void doTestBackward(Field<T> field,
391 final double epsilonLast,
392 final double epsilonMaxValue,
393 final double epsilonMaxTime,
394 final String name)
395 throws DimensionMismatchException, NumberIsTooSmallException,
396 MaxCountExceededException, NoBracketingException {
397
398 TestFieldProblem5<T> pb = new TestFieldProblem5<>(field);
399 double minStep = 0;
400 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).abs().getReal();
401 double scalAbsoluteTolerance = 1.0e-8;
402 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
403
404 EmbeddedRungeKuttaFieldIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
405 scalAbsoluteTolerance,
406 scalRelativeTolerance);
407 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<>(pb, integ);
408 integ.addStepHandler(handler);
409 integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
410
411 Assert.assertEquals(0, handler.getLastError().getReal(), epsilonLast);
412 Assert.assertEquals(0, handler.getMaximalValueError().getReal(), epsilonMaxValue);
413 Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), epsilonMaxTime);
414 Assert.assertEquals(name, integ.getName());
415 }
416
417 @Test
418 public abstract void testKepler();
419
420 protected <T extends RealFieldElement<T>> void doTestKepler(Field<T> field, double epsilon) {
421
422 final TestFieldProblem3<T> pb = new TestFieldProblem3<>(field, field.getZero().add(0.9));
423 double minStep = 0;
424 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
425 double[] vecAbsoluteTolerance = { 1.0e-8, 1.0e-8, 1.0e-10, 1.0e-10 };
426 double[] vecRelativeTolerance = { 1.0e-10, 1.0e-10, 1.0e-8, 1.0e-8 };
427
428 FirstOrderFieldIntegrator<T> integ = createIntegrator(field, minStep, maxStep,
429 vecAbsoluteTolerance, vecRelativeTolerance);
430 integ.addStepHandler(new KeplerHandler<>(pb, epsilon));
431 integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
432 }
433
434 private static final class KeplerHandler<T extends RealFieldElement<T>> implements FieldStepHandler<T> {
435 private T maxError;
436 private final TestFieldProblem3<T> pb;
437 private final double epsilon;
438 KeplerHandler(TestFieldProblem3<T> pb, double epsilon) {
439 this.pb = pb;
440 this.epsilon = epsilon;
441 maxError = pb.getField().getZero();
442 }
443 @Override
444 public void init(FieldODEStateAndDerivative<T> state0, T t) {
445 maxError = pb.getField().getZero();
446 }
447 @Override
448 public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast)
449 throws MaxCountExceededException {
450
451 FieldODEStateAndDerivative<T> current = interpolator.getCurrentState();
452 T[] theoreticalY = pb.computeTheoreticalState(current.getTime());
453 T dx = current.getState()[0].subtract(theoreticalY[0]);
454 T dy = current.getState()[1].subtract(theoreticalY[1]);
455 T error = dx.multiply(dx).add(dy.multiply(dy));
456 if (error.subtract(maxError).getReal() > 0) {
457 maxError = error;
458 }
459 if (isLast) {
460 Assert.assertEquals(0.0, maxError.getReal(), epsilon);
461 }
462 }
463 }
464
465 @Test
466 public abstract void testPartialDerivatives();
467
468 protected <T extends RealFieldElement<T>> void doTestPartialDerivatives(final double epsilonY,
469 final double[] epsilonPartials) {
470
471
472 final int parameters = 5;
473 final int order = 1;
474 final int parOmega = 0;
475 final int parTO = 1;
476 final int parY00 = 2;
477 final int parY01 = 3;
478 final int parT = 4;
479
480 DerivativeStructure omega = new DerivativeStructure(parameters, order, parOmega, 1.3);
481 DerivativeStructure t0 = new DerivativeStructure(parameters, order, parTO, 1.3);
482 DerivativeStructure[] y0 = new DerivativeStructure[] {
483 new DerivativeStructure(parameters, order, parY00, 3.0),
484 new DerivativeStructure(parameters, order, parY01, 4.0)
485 };
486 DerivativeStructure t = new DerivativeStructure(parameters, order, parT, 6.0);
487 SinCos sinCos = new SinCos(omega);
488
489 EmbeddedRungeKuttaFieldIntegrator<DerivativeStructure> integrator =
490 createIntegrator(omega.getField(),
491 t.subtract(t0).multiply(0.001).getReal(), t.subtract(t0).getReal(),
492 1.0e-12, 1.0e-12);
493 FieldODEStateAndDerivative<DerivativeStructure> result =
494 integrator.integrate(new FieldExpandableODE<>(sinCos),
495 new FieldODEState<>(t0, y0),
496 t);
497
498
499 for (int i = 0; i < sinCos.getDimension(); ++i) {
500 Assert.assertEquals(sinCos.theoreticalY(t.getReal())[i], result.getState()[i].getValue(), epsilonY);
501 }
502
503
504 final double[][] derivatives = sinCos.getDerivatives(t.getReal());
505 for (int i = 0; i < sinCos.getDimension(); ++i) {
506 for (int parameter = 0; parameter < parameters; ++parameter) {
507 Assert.assertEquals(derivatives[i][parameter], dYdP(result.getState()[i], parameter), epsilonPartials[parameter]);
508 }
509 }
510 }
511
512 private double dYdP(final DerivativeStructure y, final int parameter) {
513 int[] orders = new int[y.getFreeParameters()];
514 orders[parameter] = 1;
515 return y.getPartialDerivative(orders);
516 }
517
518 private static final class SinCos implements FirstOrderFieldDifferentialEquations<DerivativeStructure> {
519
520 private final DerivativeStructure omega;
521 private DerivativeStructure r;
522 private DerivativeStructure alpha;
523
524 private double dRdY00;
525 private double dRdY01;
526 private double dAlphadOmega;
527 private double dAlphadT0;
528 private double dAlphadY00;
529 private double dAlphadY01;
530
531 protected SinCos(final DerivativeStructure omega) {
532 this.omega = omega;
533 }
534
535 @Override
536 public int getDimension() {
537 return 2;
538 }
539
540 @Override
541 public void init(final DerivativeStructure t0, final DerivativeStructure[] y0,
542 final DerivativeStructure finalTime) {
543
544
545
546 final DerivativeStructure r2 = y0[0].multiply(y0[0]).add(y0[1].multiply(y0[1]));
547
548 this.r = r2.sqrt();
549 this.dRdY00 = y0[0].divide(r).getReal();
550 this.dRdY01 = y0[1].divide(r).getReal();
551
552 this.alpha = y0[0].atan2(y0[1]).subtract(t0.multiply(omega));
553 this.dAlphadOmega = -t0.getReal();
554 this.dAlphadT0 = -omega.getReal();
555 this.dAlphadY00 = y0[1].divide(r2).getReal();
556 this.dAlphadY01 = y0[0].negate().divide(r2).getReal();
557 }
558
559 @Override
560 public DerivativeStructure[] computeDerivatives(final DerivativeStructure t, final DerivativeStructure[] y) {
561 return new DerivativeStructure[] {
562 omega.multiply(y[1]),
563 omega.multiply(y[0]).negate()
564 };
565 }
566
567 public double[] theoreticalY(final double t) {
568 final double theta = omega.getReal() * t + alpha.getReal();
569 return new double[] {
570 r.getReal() * JdkMath.sin(theta), r.getReal() * JdkMath.cos(theta)
571 };
572 }
573
574 public double[][] getDerivatives(final double t) {
575
576
577 final double theta = omega.getReal() * t + alpha.getReal();
578 final double sin = JdkMath.sin(theta);
579 final double cos = JdkMath.cos(theta);
580 final double y0 = r.getReal() * sin;
581 final double y1 = r.getReal() * cos;
582
583
584 final double dY0dOmega = y1 * (t + dAlphadOmega);
585 final double dY0dT0 = y1 * dAlphadT0;
586 final double dY0dY00 = dRdY00 * sin + y1 * dAlphadY00;
587 final double dY0dY01 = dRdY01 * sin + y1 * dAlphadY01;
588 final double dY0dT = y1 * omega.getReal();
589
590
591 final double dY1dOmega = - y0 * (t + dAlphadOmega);
592 final double dY1dT0 = - y0 * dAlphadT0;
593 final double dY1dY00 = dRdY00 * cos - y0 * dAlphadY00;
594 final double dY1dY01 = dRdY01 * cos - y0 * dAlphadY01;
595 final double dY1dT = - y0 * omega.getReal();
596
597 return new double[][] {
598 { dY0dOmega, dY0dT0, dY0dY00, dY0dY01, dY0dT },
599 { dY1dOmega, dY1dT0, dY1dY00, dY1dY01, dY1dT }
600 };
601 }
602 }
603 }