1 /*
2  * Licensed to the Apache Software Foundation (ASF) under one or more
3  * contributor license agreements.  See the NOTICE file distributed with
4  * this work for additional information regarding copyright ownership.
5  * The ASF licenses this file to You under the Apache License, Version 2.0
6  * (the "License"); you may not use this file except in compliance with
7  * the License.  You may obtain a copy of the License at
8  *
9  *      http://www.apache.org/licenses/LICENSE-2.0
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  */
17 
18 package org.apache.commons.math3.ode;
19 
20 import org.apache.commons.math3.exception.DimensionMismatchException;
21 import org.apache.commons.math3.exception.MaxCountExceededException;
22 import org.apache.commons.math3.exception.NoBracketingException;
23 import org.apache.commons.math3.exception.NumberIsTooSmallException;
24 import org.apache.commons.math3.ode.JacobianMatrices.MismatchedEquations;
25 import org.apache.commons.math3.ode.nonstiff.DormandPrince54Integrator;
26 import org.apache.commons.math3.stat.descriptive.SummaryStatistics;
27 import org.apache.commons.math3.util.FastMath;
28 import org.junit.Assert;
29 import org.junit.Test;
30 
31 public class JacobianMatricesTest {
32 
33     @Test
testLowAccuracyExternalDifferentiation()34     public void testLowAccuracyExternalDifferentiation()
35         throws NumberIsTooSmallException, DimensionMismatchException,
36                MaxCountExceededException, NoBracketingException {
37         // this test does not really test JacobianMatrices,
38         // it only shows that WITHOUT this class, attempting to recover
39         // the jacobians from external differentiation on simple integration
40         // results with low accuracy gives very poor results. In fact,
41         // the curves dy/dp = g(b) when b varies from 2.88 to 3.08 are
42         // essentially noise.
43         // This test is taken from Hairer, Norsett and Wanner book
44         // Solving Ordinary Differential Equations I (Nonstiff problems),
45         // the curves dy/dp = g(b) are in figure 6.5
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
testHighAccuracyExternalDifferentiation()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 //            Assert.assertEquals(5000, integ.getMaxEvaluations());
141 //            Assert.assertTrue(integ.getEvaluations() > 1500);
142 //            Assert.assertTrue(integ.getEvaluations() < 2100);
143 //            Assert.assertEquals(4 * integ.getEvaluations(), integ.getEvaluations());
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 //            Assert.assertEquals(5000, integ.getMaxEvaluations());
182 //            Assert.assertTrue(integ.getEvaluations() > 350);
183 //            Assert.assertTrue(integ.getEvaluations() < 510);
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 * FastMath.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 * FastMath.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 class Brusselator extends AbstractParameterizable
311         implements MainStateJacobianProvider, ParameterJacobianProvider {
312 
313         public static final String B = "b";
314 
315         private double b;
316 
317         public Brusselator(double b) {
318             super(B);
319             this.b = b;
320         }
321 
322         public int getDimension() {
323             return 2;
324         }
325 
326         public void computeDerivatives(double t, double[] y, double[] yDot) {
327             double prod = y[0] * y[0] * y[1];
328             yDot[0] = 1 + prod - (b + 1) * y[0];
329             yDot[1] = b * y[0] - prod;
330         }
331 
332         public void computeMainStateJacobian(double t, double[] y, double[] yDot,
333                                              double[][] dFdY) {
334             double p = 2 * y[0] * y[1];
335             double y02 = y[0] * y[0];
336             dFdY[0][0] = p - (1 + b);
337             dFdY[0][1] = y02;
338             dFdY[1][0] = b - p;
339             dFdY[1][1] = -y02;
340         }
341 
342         public void computeParameterJacobian(double t, double[] y, double[] yDot,
343                                              String paramName, double[] dFdP) {
344             if (isSupported(paramName)) {
345                 dFdP[0] = -y[0];
346                 dFdP[1] = y[0];
347             } else {
348                 dFdP[0] = 0;
349                 dFdP[1] = 0;
350             }
351         }
352 
353         public double dYdP0() {
354             return -1088.232716447743 + (1050.775747149553 + (-339.012934631828 + 36.52917025056327 * b) * b) * b;
355         }
356 
357         public double dYdP1() {
358             return 1502.824469929139 + (-1438.6974831849952 + (460.959476642384 - 49.43847385647082 * b) * b) * b;
359         }
360 
361     }
362 
363     private static class ParamBrusselator extends AbstractParameterizable
364         implements FirstOrderDifferentialEquations, ParameterizedODE {
365 
366         public static final String B = "b";
367 
368         private double b;
369 
370         public ParamBrusselator(double b) {
371             super(B);
372             this.b = b;
373         }
374 
375         public int getDimension() {
376             return 2;
377         }
378 
379         /** {@inheritDoc} */
380         public double getParameter(final String name)
381             throws UnknownParameterException {
382             complainIfNotSupported(name);
383             return b;
384         }
385 
386         /** {@inheritDoc} */
387         public void setParameter(final String name, final double value)
388             throws UnknownParameterException {
389             complainIfNotSupported(name);
390             b = value;
391         }
392 
393         public void computeDerivatives(double t, double[] y, double[] yDot) {
394             double prod = y[0] * y[0] * y[1];
395             yDot[0] = 1 + prod - (b + 1) * y[0];
396             yDot[1] = b * y[0] - prod;
397         }
398 
399         public double dYdP0() {
400             return -1088.232716447743 + (1050.775747149553 + (-339.012934631828 + 36.52917025056327 * b) * b) * b;
401         }
402 
403         public double dYdP1() {
404             return 1502.824469929139 + (-1438.6974831849952 + (460.959476642384 - 49.43847385647082 * b) * b) * b;
405         }
406 
407     }
408 
409     /** ODE representing a point moving on a circle with provided center and angular rate. */
410     private static class Circle extends AbstractParameterizable
411         implements MainStateJacobianProvider, ParameterJacobianProvider {
412 
413         public static final String CX = "cx";
414         public static final String CY = "cy";
415         public static final String OMEGA = "omega";
416 
417         private final double[] y0;
418         private double cx;
419         private double cy;
420         private double omega;
421 
422         public Circle(double[] y0, double cx, double cy, double omega) {
423             super(CX,CY,OMEGA);
424             this.y0    = y0.clone();
425             this.cx    = cx;
426             this.cy    = cy;
427             this.omega = omega;
428         }
429 
430         public int getDimension() {
431             return 2;
432         }
433 
434         public void computeDerivatives(double t, double[] y, double[] yDot) {
435             yDot[0] = omega * (cy - y[1]);
436             yDot[1] = omega * (y[0] - cx);
437         }
438 
439         public void computeMainStateJacobian(double t, double[] y,
440                                              double[] yDot, double[][] dFdY) {
441             dFdY[0][0] = 0;
442             dFdY[0][1] = -omega;
443             dFdY[1][0] = omega;
444             dFdY[1][1] = 0;
445         }
446 
447         public void computeParameterJacobian(double t, double[] y, double[] yDot,
448                                              String paramName, double[] dFdP)
449             throws UnknownParameterException {
450             complainIfNotSupported(paramName);
451             if (paramName.equals(CX)) {
452                 dFdP[0] = 0;
453                 dFdP[1] = -omega;
454             } else if (paramName.equals(CY)) {
455                 dFdP[0] = omega;
456                 dFdP[1] = 0;
457             }  else {
458                 dFdP[0] = cy - y[1];
459                 dFdP[1] = y[0] - cx;
460             }
461         }
462 
463         public double[] exactY(double t) {
464             double cos = FastMath.cos(omega * t);
465             double sin = FastMath.sin(omega * t);
466             double dx0 = y0[0] - cx;
467             double dy0 = y0[1] - cy;
468             return new double[] {
469                 cx + cos * dx0 - sin * dy0,
470                 cy + sin * dx0 + cos * dy0
471             };
472         }
473 
474         public double[][] exactDyDy0(double t) {
475             double cos = FastMath.cos(omega * t);
476             double sin = FastMath.sin(omega * t);
477             return new double[][] {
478                 { cos, -sin },
479                 { sin,  cos }
480             };
481         }
482 
483         public double[] exactDyDcx(double t) {
484             double cos = FastMath.cos(omega * t);
485             double sin = FastMath.sin(omega * t);
486             return new double[] {1 - cos, -sin};
487         }
488 
489         public double[] exactDyDcy(double t) {
490             double cos = FastMath.cos(omega * t);
491             double sin = FastMath.sin(omega * t);
492             return new double[] {sin, 1 - cos};
493         }
494 
495         public double[] exactDyDom(double t) {
496             double cos = FastMath.cos(omega * t);
497             double sin = FastMath.sin(omega * t);
498             double dx0 = y0[0] - cx;
499             double dy0 = y0[1] - cy;
500             return new double[] { -t * (sin * dx0 + cos * dy0) , t * (cos * dx0 - sin * dy0) };
501         }
502 
503     }
504 
505     /** ODE representing a point moving on a circle with provided center and angular rate. */
506     private static class ParameterizedCircle extends AbstractParameterizable
507         implements FirstOrderDifferentialEquations, ParameterizedODE {
508 
509         public static final String CX = "cx";
510         public static final String CY = "cy";
511         public static final String OMEGA = "omega";
512 
513         private final double[] y0;
514         private double cx;
515         private double cy;
516         private double omega;
517 
518         public ParameterizedCircle(double[] y0, double cx, double cy, double omega) {
519             super(CX,CY,OMEGA);
520             this.y0    = y0.clone();
521             this.cx    = cx;
522             this.cy    = cy;
523             this.omega = omega;
524         }
525 
526         public int getDimension() {
527             return 2;
528         }
529 
530         public void computeDerivatives(double t, double[] y, double[] yDot) {
531             yDot[0] = omega * (cy - y[1]);
532             yDot[1] = omega * (y[0] - cx);
533         }
534 
535         public double getParameter(final String name)
536             throws UnknownParameterException {
537             if (name.equals(CX)) {
538                 return cx;
539             } else if (name.equals(CY)) {
540                     return cy;
541             } else if (name.equals(OMEGA)) {
542                 return omega;
543             } else {
544                 throw new UnknownParameterException(name);
545             }
546         }
547 
548         public void setParameter(final String name, final double value)
549             throws UnknownParameterException {
550             if (name.equals(CX)) {
551                 cx = value;
552             } else if (name.equals(CY)) {
553                 cy = value;
554             } else if (name.equals(OMEGA)) {
555                 omega = value;
556             } else {
557                 throw new UnknownParameterException(name);
558             }
559         }
560 
561         public double[] exactY(double t) {
562             double cos = FastMath.cos(omega * t);
563             double sin = FastMath.sin(omega * t);
564             double dx0 = y0[0] - cx;
565             double dy0 = y0[1] - cy;
566             return new double[] {
567                 cx + cos * dx0 - sin * dy0,
568                 cy + sin * dx0 + cos * dy0
569             };
570         }
571 
572         public double[][] exactDyDy0(double t) {
573             double cos = FastMath.cos(omega * t);
574             double sin = FastMath.sin(omega * t);
575             return new double[][] {
576                 { cos, -sin },
577                 { sin,  cos }
578             };
579         }
580 
581         public double[] exactDyDcx(double t) {
582             double cos = FastMath.cos(omega * t);
583             double sin = FastMath.sin(omega * t);
584             return new double[] {1 - cos, -sin};
585         }
586 
587         public double[] exactDyDcy(double t) {
588             double cos = FastMath.cos(omega * t);
589             double sin = FastMath.sin(omega * t);
590             return new double[] {sin, 1 - cos};
591         }
592 
593         public double[] exactDyDom(double t) {
594             double cos = FastMath.cos(omega * t);
595             double sin = FastMath.sin(omega * t);
596             double dx0 = y0[0] - cx;
597             double dy0 = y0[1] - cy;
598             return new double[] { -t * (sin * dx0 + cos * dy0) , t * (cos * dx0 - sin * dy0) };
599         }
600 
601     }
602 
603 }
604