1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 package org.hipparchus.migration.ode;
24
25 import org.hipparchus.exception.MathIllegalArgumentException;
26 import org.hipparchus.exception.MathIllegalStateException;
27 import org.hipparchus.ode.AbstractIntegrator;
28 import org.hipparchus.ode.AbstractParameterizable;
29 import org.hipparchus.ode.ExpandableODE;
30 import org.hipparchus.ode.LocalizedODEFormats;
31 import org.hipparchus.ode.NamedParameterJacobianProvider;
32 import org.hipparchus.ode.ODEIntegrator;
33 import org.hipparchus.ode.ODEState;
34 import org.hipparchus.ode.ODEStateAndDerivative;
35 import org.hipparchus.ode.OrdinaryDifferentialEquation;
36 import org.hipparchus.ode.ParametersController;
37 import org.hipparchus.ode.nonstiff.DormandPrince54Integrator;
38 import org.hipparchus.stat.descriptive.StreamingStatistics;
39 import org.hipparchus.util.FastMath;
40 import org.junit.Assert;
41 import org.junit.Test;
42
43 @Deprecated
44 public class JacobianMatricesTest {
45
46 @Test
47 public void testLowAccuracyExternalDifferentiation()
48 throws MathIllegalArgumentException, MathIllegalStateException {
49
50
51
52
53
54
55
56
57
58 ODEIntegrator integ =
59 new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 });
60 double hP = 1.0e-12;
61 StreamingStatistics residualsP0 = new StreamingStatistics();
62 StreamingStatistics residualsP1 = new StreamingStatistics();
63 for (double b = 2.88; b < 3.08; b += 0.001) {
64 Brusselator brusselator = new Brusselator(b);
65 double[] y = { 1.3, b };
66 y = integ.integrate(brusselator, new ODEState(0, y), 20.0).getPrimaryState();
67 double[] yP = { 1.3, b + hP };
68 yP = integ.integrate(brusselator, new ODEState(0, yP), 20.0).getPrimaryState();
69 residualsP0.addValue((yP[0] - y[0]) / hP - brusselator.dYdP0());
70 residualsP1.addValue((yP[1] - y[1]) / hP - brusselator.dYdP1());
71 }
72 Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) > 500);
73 Assert.assertTrue(residualsP0.getStandardDeviation() > 30);
74 Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) > 700);
75 Assert.assertTrue(residualsP1.getStandardDeviation() > 40);
76 }
77
78 @Test
79 public void testHighAccuracyExternalDifferentiation()
80 throws MathIllegalArgumentException, MathIllegalStateException {
81 ODEIntegrator integ =
82 new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 1.0e-10, 1.0e-10 });
83 double hP = 1.0e-12;
84 StreamingStatistics residualsP0 = new StreamingStatistics();
85 StreamingStatistics residualsP1 = new StreamingStatistics();
86 for (double b = 2.88; b < 3.08; b += 0.001) {
87 ParamBrusselator brusselator = new ParamBrusselator(b);
88 double[] y = { 1.3, b };
89 y = integ.integrate(brusselator, new ODEState(0, y), 20.0).getPrimaryState();
90 double[] yP = { 1.3, b + hP };
91 brusselator.setParameter("b", b + hP);
92 yP = integ.integrate(brusselator, new ODEState(0, yP), 20.0).getPrimaryState();
93 residualsP0.addValue((yP[0] - y[0]) / hP - brusselator.dYdP0());
94 residualsP1.addValue((yP[1] - y[1]) / hP - brusselator.dYdP1());
95 }
96 Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) > 0.02);
97 Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.03);
98 Assert.assertTrue(residualsP0.getStandardDeviation() > 0.003);
99 Assert.assertTrue(residualsP0.getStandardDeviation() < 0.004);
100 Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) > 0.04);
101 Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05);
102 Assert.assertTrue(residualsP1.getStandardDeviation() > 0.007);
103 Assert.assertTrue(residualsP1.getStandardDeviation() < 0.008);
104 }
105
106 @Test
107 public void testWrongParameterName() {
108 final String name = "an-unknown-parameter";
109 try {
110 ParamBrusselator brusselator = new ParamBrusselator(2.9);
111 brusselator.setParameter(name, 3.0);
112 Assert.fail("an exception should have been thrown");
113 } catch (MathIllegalArgumentException upe) {
114 Assert.assertEquals(LocalizedODEFormats.UNKNOWN_PARAMETER, upe.getSpecifier());
115 Assert.assertEquals(name, upe.getParts()[0]);
116 }
117 }
118
119 @Test
120 public void testInternalDifferentiation()
121 throws MathIllegalArgumentException, MathIllegalStateException {
122 AbstractIntegrator integ =
123 new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 });
124 double hP = 1.0e-12;
125 double hY = 1.0e-12;
126 StreamingStatistics residualsP0 = new StreamingStatistics();
127 StreamingStatistics residualsP1 = new StreamingStatistics();
128 for (double b = 2.88; b < 3.08; b += 0.001) {
129 ParamBrusselator brusselator = new ParamBrusselator(b);
130 brusselator.setParameter(ParamBrusselator.B, b);
131 double[] z = { 1.3, b };
132
133 JacobianMatrices jacob = new JacobianMatrices(brusselator, new double[] { hY, hY }, ParamBrusselator.B);
134 jacob.setParametersController(brusselator);
135 jacob.setParameterStep(ParamBrusselator.B, hP);
136 jacob.setInitialParameterJacobian(ParamBrusselator.B, new double[] { 0.0, 1.0 });
137
138 ExpandableODE efode = new ExpandableODE(brusselator);
139 jacob.registerVariationalEquations(efode);
140
141 integ.setMaxEvaluations(5000);
142 final ODEState initialState = jacob.setUpInitialState(new ODEState(0, z));
143 final ODEStateAndDerivative finalState = integ.integrate(efode, initialState, 20.0);
144 final double[] dZdP = jacob.extractParameterJacobian(finalState, ParamBrusselator.B);
145
146
147
148
149 residualsP0.addValue(dZdP[0] - brusselator.dYdP0());
150 residualsP1.addValue(dZdP[1] - brusselator.dYdP1());
151 }
152 Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.02);
153 Assert.assertTrue(residualsP0.getStandardDeviation() < 0.003);
154 Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05);
155 Assert.assertTrue(residualsP1.getStandardDeviation() < 0.01);
156 }
157
158 @Test
159 public void testAnalyticalDifferentiation()
160 throws MathIllegalArgumentException, MathIllegalStateException {
161 AbstractIntegrator integ =
162 new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-4, 1.0e-4 }, new double[] { 1.0e-4, 1.0e-4 });
163 StreamingStatistics residualsP0 = new StreamingStatistics();
164 StreamingStatistics residualsP1 = new StreamingStatistics();
165 for (double b = 2.88; b < 3.08; b += 0.001) {
166 Brusselator brusselator = new Brusselator(b);
167 double[] z = { 1.3, b };
168
169 JacobianMatrices jacob = new JacobianMatrices(brusselator, Brusselator.B);
170 jacob.addParameterJacobianProvider(brusselator);
171 jacob.setInitialParameterJacobian(Brusselator.B, new double[] { 0.0, 1.0 });
172
173 ExpandableODE efode = new ExpandableODE(brusselator);
174 jacob.registerVariationalEquations(efode);
175
176 integ.setMaxEvaluations(5000);
177 final ODEState initialState = jacob.setUpInitialState(new ODEState(0, z));
178 final ODEStateAndDerivative finalState = integ.integrate(efode, initialState, 20.0);
179 final double[] dZdP = jacob.extractParameterJacobian(finalState, Brusselator.B);
180
181
182
183 residualsP0.addValue(dZdP[0] - brusselator.dYdP0());
184 residualsP1.addValue(dZdP[1] - brusselator.dYdP1());
185 }
186 Assert.assertTrue((residualsP0.getMax() - residualsP0.getMin()) < 0.014);
187 Assert.assertTrue(residualsP0.getStandardDeviation() < 0.003);
188 Assert.assertTrue((residualsP1.getMax() - residualsP1.getMin()) < 0.05);
189 Assert.assertTrue(residualsP1.getStandardDeviation() < 0.01);
190 }
191
192 @Test
193 public void testFinalResult()
194 throws MathIllegalArgumentException, MathIllegalStateException {
195
196 AbstractIntegrator integ =
197 new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 1.0e-10, 1.0e-10 });
198 double[] y = new double[] { 0.0, 1.0 };
199 Circle circle = new Circle(y, 1.0, 1.0, 0.1);
200
201 JacobianMatrices jacob = new JacobianMatrices(circle, Circle.CX, Circle.CY, Circle.OMEGA);
202 jacob.addParameterJacobianProvider(circle);
203 jacob.setInitialMainStateJacobian(circle.exactDyDy0(0));
204 jacob.setInitialParameterJacobian(Circle.CX, circle.exactDyDcx(0));
205 jacob.setInitialParameterJacobian(Circle.CY, circle.exactDyDcy(0));
206 jacob.setInitialParameterJacobian(Circle.OMEGA, circle.exactDyDom(0));
207
208 ExpandableODE efode = new ExpandableODE(circle);
209 jacob.registerVariationalEquations(efode);
210
211 integ.setMaxEvaluations(5000);
212
213 double t = 18 * FastMath.PI;
214 final ODEState initialState = jacob.setUpInitialState(new ODEState(0, y));
215 final ODEStateAndDerivative finalState = integ.integrate(efode, initialState, t);
216 y = finalState.getPrimaryState();
217 for (int i = 0; i < y.length; ++i) {
218 Assert.assertEquals(circle.exactY(t)[i], y[i], 1.0e-9);
219 }
220
221 double[][] dydy0 = jacob.extractMainSetJacobian(finalState);
222 for (int i = 0; i < dydy0.length; ++i) {
223 for (int j = 0; j < dydy0[i].length; ++j) {
224 Assert.assertEquals(circle.exactDyDy0(t)[i][j], dydy0[i][j], 1.0e-9);
225 }
226 }
227 double[] dydcx = jacob.extractParameterJacobian(finalState, Circle.CX);
228 for (int i = 0; i < dydcx.length; ++i) {
229 Assert.assertEquals(circle.exactDyDcx(t)[i], dydcx[i], 1.0e-7);
230 }
231 double[] dydcy = jacob.extractParameterJacobian(finalState, Circle.CY);
232 for (int i = 0; i < dydcy.length; ++i) {
233 Assert.assertEquals(circle.exactDyDcy(t)[i], dydcy[i], 1.0e-7);
234 }
235 double[] dydom = jacob.extractParameterJacobian(finalState, Circle.OMEGA);
236 for (int i = 0; i < dydom.length; ++i) {
237 Assert.assertEquals(circle.exactDyDom(t)[i], dydom[i], 1.0e-7);
238 }
239 }
240
241 @Test
242 public void testParameterizable()
243 throws MathIllegalArgumentException, MathIllegalStateException {
244
245 AbstractIntegrator integ =
246 new DormandPrince54Integrator(1.0e-8, 100.0, new double[] { 1.0e-10, 1.0e-10 }, new double[] { 1.0e-10, 1.0e-10 });
247 double[] y = new double[] { 0.0, 1.0 };
248 ParameterizedCircle pcircle = new ParameterizedCircle(y, 1.0, 1.0, 0.1);
249
250 double hP = 1.0e-12;
251 double hY = 1.0e-12;
252
253 JacobianMatrices jacob = new JacobianMatrices(pcircle, new double[] { hY, hY },
254 ParameterizedCircle.CX, ParameterizedCircle.CY,
255 ParameterizedCircle.OMEGA);
256 jacob.setParametersController(pcircle);
257 jacob.setParameterStep(ParameterizedCircle.CX, hP);
258 jacob.setParameterStep(ParameterizedCircle.CY, hP);
259 jacob.setParameterStep(ParameterizedCircle.OMEGA, hP);
260 jacob.setInitialMainStateJacobian(pcircle.exactDyDy0(0));
261 jacob.setInitialParameterJacobian(ParameterizedCircle.CX, pcircle.exactDyDcx(0));
262 jacob.setInitialParameterJacobian(ParameterizedCircle.CY, pcircle.exactDyDcy(0));
263 jacob.setInitialParameterJacobian(ParameterizedCircle.OMEGA, pcircle.exactDyDom(0));
264
265 ExpandableODE efode = new ExpandableODE(pcircle);
266 jacob.registerVariationalEquations(efode);
267
268 integ.setMaxEvaluations(50000);
269
270 double t = 18 * FastMath.PI;
271 final ODEState initialState = jacob.setUpInitialState(new ODEState(0, y));
272 final ODEStateAndDerivative finalState = integ.integrate(efode, initialState, t);
273 y = finalState.getPrimaryState();
274 for (int i = 0; i < y.length; ++i) {
275 Assert.assertEquals(pcircle.exactY(t)[i], y[i], 1.0e-9);
276 }
277
278 double[][] dydy0 = jacob.extractMainSetJacobian(finalState);
279 for (int i = 0; i < dydy0.length; ++i) {
280 for (int j = 0; j < dydy0[i].length; ++j) {
281 Assert.assertEquals(pcircle.exactDyDy0(t)[i][j], dydy0[i][j], 5.0e-4);
282 }
283 }
284
285 double[] dydp0 = jacob.extractParameterJacobian(finalState, ParameterizedCircle.CX);
286 for (int i = 0; i < dydp0.length; ++i) {
287 Assert.assertEquals(pcircle.exactDyDcx(t)[i], dydp0[i], 5.0e-4);
288 }
289
290 double[] dydp1 = jacob.extractParameterJacobian(finalState, ParameterizedCircle.OMEGA);
291 for (int i = 0; i < dydp1.length; ++i) {
292 Assert.assertEquals(pcircle.exactDyDom(t)[i], dydp1[i], 1.0e-2);
293 }
294 }
295
296 private static class Brusselator extends AbstractParameterizable
297 implements MainStateJacobianProvider, NamedParameterJacobianProvider {
298
299 public static final String B = "b";
300
301 private double b;
302
303 public Brusselator(double b) {
304 super(B);
305 this.b = b;
306 }
307
308 @Override
309 public int getDimension() {
310 return 2;
311 }
312
313 @Override
314 public double[] computeDerivatives(double t, double[] y) {
315 double prod = y[0] * y[0] * y[1];
316 return new double[] {
317 1 + prod - (b + 1) * y[0],
318 b * y[0] - prod
319 };
320 }
321
322 @Override
323 public double[][] computeMainStateJacobian(double t, double[] y, double[] yDot) {
324 double p = 2 * y[0] * y[1];
325 double y02 = y[0] * y[0];
326 return new double[][] {
327 { p - (1 + b), y02 },
328 { b - p, -y02}
329 };
330 }
331
332 @Override
333 public double[] computeParameterJacobian(double t, double[] y, double[] yDot,
334 String paramName) {
335 if (isSupported(paramName)) {
336 return new double[] { -y[0], y[0] };
337 } else {
338 return new double[] { 0, 0 };
339 }
340 }
341
342 public double dYdP0() {
343 return -1088.232716447743 + (1050.775747149553 + (-339.012934631828 + 36.52917025056327 * b) * b) * b;
344 }
345
346 public double dYdP1() {
347 return 1502.824469929139 + (-1438.6974831849952 + (460.959476642384 - 49.43847385647082 * b) * b) * b;
348 }
349
350 }
351
352 private static class ParamBrusselator extends AbstractParameterizable
353 implements OrdinaryDifferentialEquation, ParametersController {
354
355 public static final String B = "b";
356
357 private double b;
358
359 public ParamBrusselator(double b) {
360 super(B);
361 this.b = b;
362 }
363
364 @Override
365 public int getDimension() {
366 return 2;
367 }
368
369
370 @Override
371 public double getParameter(final String name)
372 throws MathIllegalArgumentException {
373 complainIfNotSupported(name);
374 return b;
375 }
376
377
378 @Override
379 public void setParameter(final String name, final double value)
380 throws MathIllegalArgumentException {
381 complainIfNotSupported(name);
382 b = value;
383 }
384
385 @Override
386 public double[] computeDerivatives(double t, double[] y) {
387 double prod = y[0] * y[0] * y[1];
388 return new double[] {
389 1 + prod - (b + 1) * y[0],
390 b * y[0] - prod
391 };
392 }
393
394 public double dYdP0() {
395 return -1088.232716447743 + (1050.775747149553 + (-339.012934631828 + 36.52917025056327 * b) * b) * b;
396 }
397
398 public double dYdP1() {
399 return 1502.824469929139 + (-1438.6974831849952 + (460.959476642384 - 49.43847385647082 * b) * b) * b;
400 }
401
402 }
403
404
405 private static class Circle extends AbstractParameterizable
406 implements MainStateJacobianProvider, NamedParameterJacobianProvider {
407
408 public static final String CX = "cx";
409 public static final String CY = "cy";
410 public static final String OMEGA = "omega";
411
412 private final double[] y0;
413 private double cx;
414 private double cy;
415 private double omega;
416
417 public Circle(double[] y0, double cx, double cy, double omega) {
418 super(CX,CY,OMEGA);
419 this.y0 = y0.clone();
420 this.cx = cx;
421 this.cy = cy;
422 this.omega = omega;
423 }
424
425 @Override
426 public int getDimension() {
427 return 2;
428 }
429
430 @Override
431 public double[] computeDerivatives(double t, double[] y) {
432 return new double[] {
433 omega * (cy - y[1]),
434 omega * (y[0] - cx)
435 };
436 }
437
438 @Override
439 public double[][] computeMainStateJacobian(double t, double[] y, double[] yDot) {
440 return new double[][] {
441 { 0, -omega },
442 { omega, 0 }
443 };
444 }
445
446 @Override
447 public double[] computeParameterJacobian(double t, double[] y, double[] yDot, String paramName)
448 throws MathIllegalArgumentException {
449 complainIfNotSupported(paramName);
450 if (paramName.equals(CX)) {
451 return new double[] { 0, -omega };
452 } else if (paramName.equals(CY)) {
453 return new double[] { omega, 0 };
454 } else {
455 return new double[] { cy - y[1], y[0] - cx };
456 }
457 }
458
459 public double[] exactY(double t) {
460 double cos = FastMath.cos(omega * t);
461 double sin = FastMath.sin(omega * t);
462 double dx0 = y0[0] - cx;
463 double dy0 = y0[1] - cy;
464 return new double[] {
465 cx + cos * dx0 - sin * dy0,
466 cy + sin * dx0 + cos * dy0
467 };
468 }
469
470 public double[][] exactDyDy0(double t) {
471 double cos = FastMath.cos(omega * t);
472 double sin = FastMath.sin(omega * t);
473 return new double[][] {
474 { cos, -sin },
475 { sin, cos }
476 };
477 }
478
479 public double[] exactDyDcx(double t) {
480 double cos = FastMath.cos(omega * t);
481 double sin = FastMath.sin(omega * t);
482 return new double[] {1 - cos, -sin};
483 }
484
485 public double[] exactDyDcy(double t) {
486 double cos = FastMath.cos(omega * t);
487 double sin = FastMath.sin(omega * t);
488 return new double[] {sin, 1 - cos};
489 }
490
491 public double[] exactDyDom(double t) {
492 double cos = FastMath.cos(omega * t);
493 double sin = FastMath.sin(omega * t);
494 double dx0 = y0[0] - cx;
495 double dy0 = y0[1] - cy;
496 return new double[] { -t * (sin * dx0 + cos * dy0) , t * (cos * dx0 - sin * dy0) };
497 }
498
499 }
500
501
502 private static class ParameterizedCircle extends AbstractParameterizable
503 implements OrdinaryDifferentialEquation, ParametersController {
504
505 public static final String CX = "cx";
506 public static final String CY = "cy";
507 public static final String OMEGA = "omega";
508
509 private final double[] y0;
510 private double cx;
511 private double cy;
512 private double omega;
513
514 public ParameterizedCircle(double[] y0, double cx, double cy, double omega) {
515 super(CX,CY,OMEGA);
516 this.y0 = y0.clone();
517 this.cx = cx;
518 this.cy = cy;
519 this.omega = omega;
520 }
521
522 @Override
523 public int getDimension() {
524 return 2;
525 }
526
527 @Override
528 public double[] computeDerivatives(double t, double[] y) {
529 return new double[] {
530 omega * (cy - y[1]),
531 omega * (y[0] - cx)
532 };
533 }
534
535 @Override
536 public double getParameter(final String name)
537 throws MathIllegalArgumentException {
538 if (name.equals(CX)) {
539 return cx;
540 } else if (name.equals(CY)) {
541 return cy;
542 } else if (name.equals(OMEGA)) {
543 return omega;
544 } else {
545 throw new MathIllegalArgumentException(LocalizedODEFormats.UNKNOWN_PARAMETER, name);
546 }
547 }
548
549 @Override
550 public void setParameter(final String name, final double value)
551 throws MathIllegalArgumentException {
552 if (name.equals(CX)) {
553 cx = value;
554 } else if (name.equals(CY)) {
555 cy = value;
556 } else if (name.equals(OMEGA)) {
557 omega = value;
558 } else {
559 throw new MathIllegalArgumentException(LocalizedODEFormats.UNKNOWN_PARAMETER, name);
560 }
561 }
562
563 public double[] exactY(double t) {
564 double cos = FastMath.cos(omega * t);
565 double sin = FastMath.sin(omega * t);
566 double dx0 = y0[0] - cx;
567 double dy0 = y0[1] - cy;
568 return new double[] {
569 cx + cos * dx0 - sin * dy0,
570 cy + sin * dx0 + cos * dy0
571 };
572 }
573
574 public double[][] exactDyDy0(double t) {
575 double cos = FastMath.cos(omega * t);
576 double sin = FastMath.sin(omega * t);
577 return new double[][] {
578 { cos, -sin },
579 { sin, cos }
580 };
581 }
582
583 public double[] exactDyDcx(double t) {
584 double cos = FastMath.cos(omega * t);
585 double sin = FastMath.sin(omega * t);
586 return new double[] {1 - cos, -sin};
587 }
588
589 public double[] exactDyDcy(double t) {
590 double cos = FastMath.cos(omega * t);
591 double sin = FastMath.sin(omega * t);
592 return new double[] {sin, 1 - cos};
593 }
594
595 public double[] exactDyDom(double t) {
596 double cos = FastMath.cos(omega * t);
597 double sin = FastMath.sin(omega * t);
598 double dx0 = y0[0] - cx;
599 double dy0 = y0[1] - cy;
600 return new double[] { -t * (sin * dx0 + cos * dy0) , t * (cos * dx0 - sin * dy0) };
601 }
602
603 }
604
605 }