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