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