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.nonstiff; 19 20 21 import java.lang.reflect.Array; 22 23 import org.apache.commons.math3.Field; 24 import org.apache.commons.math3.RealFieldElement; 25 import org.apache.commons.math3.analysis.differentiation.DerivativeStructure; 26 import org.apache.commons.math3.exception.DimensionMismatchException; 27 import org.apache.commons.math3.exception.MaxCountExceededException; 28 import org.apache.commons.math3.exception.NoBracketingException; 29 import org.apache.commons.math3.exception.NumberIsTooSmallException; 30 import org.apache.commons.math3.ode.FieldExpandableODE; 31 import org.apache.commons.math3.ode.FirstOrderFieldDifferentialEquations; 32 import org.apache.commons.math3.ode.FieldODEState; 33 import org.apache.commons.math3.ode.FieldODEStateAndDerivative; 34 import org.apache.commons.math3.ode.TestFieldProblem1; 35 import org.apache.commons.math3.ode.TestFieldProblem2; 36 import org.apache.commons.math3.ode.TestFieldProblem3; 37 import org.apache.commons.math3.ode.TestFieldProblem4; 38 import org.apache.commons.math3.ode.TestFieldProblem5; 39 import org.apache.commons.math3.ode.TestFieldProblem6; 40 import org.apache.commons.math3.ode.TestFieldProblemAbstract; 41 import org.apache.commons.math3.ode.TestFieldProblemHandler; 42 import org.apache.commons.math3.ode.events.Action; 43 import org.apache.commons.math3.ode.events.FieldEventHandler; 44 import org.apache.commons.math3.ode.sampling.FieldStepHandler; 45 import org.apache.commons.math3.ode.sampling.FieldStepInterpolator; 46 import org.apache.commons.math3.ode.sampling.StepInterpolatorTestUtils; 47 import org.apache.commons.math3.util.FastMath; 48 import org.apache.commons.math3.util.MathArrays; 49 import org.junit.Assert; 50 import org.junit.Test; 51 52 public abstract class RungeKuttaFieldIntegratorAbstractTest { 53 54 protected abstract <T extends RealFieldElement<T>> RungeKuttaFieldIntegrator<T> createIntegrator(Field<T> field, T step)55 createIntegrator(Field<T> field, T step); 56 57 @Test testNonFieldIntegratorConsistency()58 public abstract void testNonFieldIntegratorConsistency(); 59 doTestNonFieldIntegratorConsistency(final Field<T> field)60 protected <T extends RealFieldElement<T>> void doTestNonFieldIntegratorConsistency(final Field<T> field) { 61 try { 62 63 // get the Butcher arrays from the field integrator 64 RungeKuttaFieldIntegrator<T> fieldIntegrator = createIntegrator(field, field.getZero().add(1)); 65 T[][] fieldA = fieldIntegrator.getA(); 66 T[] fieldB = fieldIntegrator.getB(); 67 T[] fieldC = fieldIntegrator.getC(); 68 69 String fieldName = fieldIntegrator.getClass().getName(); 70 String regularName = fieldName.replaceAll("Field", ""); 71 72 // get the Butcher arrays from the regular integrator 73 @SuppressWarnings("unchecked") 74 Class<RungeKuttaIntegrator> c = (Class<RungeKuttaIntegrator>) Class.forName(regularName); 75 java.lang.reflect.Field jlrFieldA = c.getDeclaredField("STATIC_A"); 76 jlrFieldA.setAccessible(true); 77 double[][] regularA = (double[][]) jlrFieldA.get(null); 78 java.lang.reflect.Field jlrFieldB = c.getDeclaredField("STATIC_B"); 79 jlrFieldB.setAccessible(true); 80 double[] regularB = (double[]) jlrFieldB.get(null); 81 java.lang.reflect.Field jlrFieldC = c.getDeclaredField("STATIC_C"); 82 jlrFieldC.setAccessible(true); 83 double[] regularC = (double[]) jlrFieldC.get(null); 84 85 Assert.assertEquals(regularA.length, fieldA.length); 86 for (int i = 0; i < regularA.length; ++i) { 87 checkArray(regularA[i], fieldA[i]); 88 } 89 checkArray(regularB, fieldB); 90 checkArray(regularC, fieldC); 91 92 } catch (ClassNotFoundException cnfe) { 93 Assert.fail(cnfe.getLocalizedMessage()); 94 } catch (IllegalAccessException iae) { 95 Assert.fail(iae.getLocalizedMessage()); 96 } catch (IllegalArgumentException iae) { 97 Assert.fail(iae.getLocalizedMessage()); 98 } catch (SecurityException se) { 99 Assert.fail(se.getLocalizedMessage()); 100 } catch (NoSuchFieldException nsfe) { 101 Assert.fail(nsfe.getLocalizedMessage()); 102 } 103 } 104 checkArray(double[] regularArray, T[] fieldArray)105 private <T extends RealFieldElement<T>> void checkArray(double[] regularArray, T[] fieldArray) { 106 Assert.assertEquals(regularArray.length, fieldArray.length); 107 for (int i = 0; i < regularArray.length; ++i) { 108 if (regularArray[i] == 0) { 109 Assert.assertTrue(0.0 == fieldArray[i].getReal()); 110 } else { 111 Assert.assertEquals(regularArray[i], fieldArray[i].getReal(), FastMath.ulp(regularArray[i])); 112 } 113 } 114 } 115 116 @Test testMissedEndEvent()117 public abstract void testMissedEndEvent(); 118 doTestMissedEndEvent(final Field<T> field, final double epsilonT, final double epsilonY)119 protected <T extends RealFieldElement<T>> void doTestMissedEndEvent(final Field<T> field, 120 final double epsilonT, final double epsilonY) 121 throws DimensionMismatchException, NumberIsTooSmallException, 122 MaxCountExceededException, NoBracketingException { 123 final T t0 = field.getZero().add(1878250320.0000029); 124 final T tEvent = field.getZero().add(1878250379.9999986); 125 final T[] k = MathArrays.buildArray(field, 3); 126 k[0] = field.getZero().add(1.0e-4); 127 k[1] = field.getZero().add(1.0e-5); 128 k[2] = field.getZero().add(1.0e-6); 129 FirstOrderFieldDifferentialEquations<T> ode = new FirstOrderFieldDifferentialEquations<T>() { 130 131 public int getDimension() { 132 return k.length; 133 } 134 135 public void init(T t0, T[] y0, T t) { 136 } 137 138 public T[] computeDerivatives(T t, T[] y) { 139 T[] yDot = MathArrays.buildArray(field, k.length); 140 for (int i = 0; i < y.length; ++i) { 141 yDot[i] = k[i].multiply(y[i]); 142 } 143 return yDot; 144 } 145 }; 146 147 RungeKuttaFieldIntegrator<T> integrator = createIntegrator(field, field.getZero().add(60.0)); 148 149 T[] y0 = MathArrays.buildArray(field, k.length); 150 for (int i = 0; i < y0.length; ++i) { 151 y0[i] = field.getOne().add(i); 152 } 153 154 FieldODEStateAndDerivative<T> result = integrator.integrate(new FieldExpandableODE<T>(ode), 155 new FieldODEState<T>(t0, y0), 156 tEvent); 157 Assert.assertEquals(tEvent.getReal(), result.getTime().getReal(), epsilonT); 158 T[] y = result.getState(); 159 for (int i = 0; i < y.length; ++i) { 160 Assert.assertEquals(y0[i].multiply(k[i].multiply(result.getTime().subtract(t0)).exp()).getReal(), 161 y[i].getReal(), 162 epsilonY); 163 } 164 165 integrator.addEventHandler(new FieldEventHandler<T>() { 166 167 public void init(FieldODEStateAndDerivative<T> state0, T t) { 168 } 169 170 public FieldODEState<T> resetState(FieldODEStateAndDerivative<T> state) { 171 return state; 172 } 173 174 public T g(FieldODEStateAndDerivative<T> state) { 175 return state.getTime().subtract(tEvent); 176 } 177 178 public Action eventOccurred(FieldODEStateAndDerivative<T> state, boolean increasing) { 179 Assert.assertEquals(tEvent.getReal(), state.getTime().getReal(), epsilonT); 180 return Action.CONTINUE; 181 } 182 }, Double.POSITIVE_INFINITY, 1.0e-20, 100); 183 result = integrator.integrate(new FieldExpandableODE<T>(ode), 184 new FieldODEState<T>(t0, y0), 185 tEvent.add(120)); 186 Assert.assertEquals(tEvent.add(120).getReal(), result.getTime().getReal(), epsilonT); 187 y = result.getState(); 188 for (int i = 0; i < y.length; ++i) { 189 Assert.assertEquals(y0[i].multiply(k[i].multiply(result.getTime().subtract(t0)).exp()).getReal(), 190 y[i].getReal(), 191 epsilonY); 192 } 193 194 } 195 196 @Test testSanityChecks()197 public abstract void testSanityChecks(); 198 doTestSanityChecks(Field<T> field)199 protected <T extends RealFieldElement<T>> void doTestSanityChecks(Field<T> field) 200 throws DimensionMismatchException, NumberIsTooSmallException, 201 MaxCountExceededException, NoBracketingException { 202 RungeKuttaFieldIntegrator<T> integrator = createIntegrator(field, field.getZero().add(0.01)); 203 try { 204 TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field); 205 integrator.integrate(new FieldExpandableODE<T>(pb), 206 new FieldODEState<T>(field.getZero(), MathArrays.buildArray(field, pb.getDimension() + 10)), 207 field.getOne()); 208 Assert.fail("an exception should have been thrown"); 209 } catch(DimensionMismatchException ie) { 210 } 211 try { 212 TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field); 213 integrator.integrate(new FieldExpandableODE<T>(pb), 214 new FieldODEState<T>(field.getZero(), MathArrays.buildArray(field, pb.getDimension())), 215 field.getZero()); 216 Assert.fail("an exception should have been thrown"); 217 } catch(NumberIsTooSmallException ie) { 218 } 219 } 220 221 @Test testDecreasingSteps()222 public abstract void testDecreasingSteps(); 223 doTestDecreasingSteps(Field<T> field, final double safetyValueFactor, final double safetyTimeFactor, final double epsilonT)224 protected <T extends RealFieldElement<T>> void doTestDecreasingSteps(Field<T> field, 225 final double safetyValueFactor, 226 final double safetyTimeFactor, 227 final double epsilonT) 228 throws DimensionMismatchException, NumberIsTooSmallException, 229 MaxCountExceededException, NoBracketingException { 230 231 @SuppressWarnings("unchecked") 232 TestFieldProblemAbstract<T>[] allProblems = 233 (TestFieldProblemAbstract<T>[]) Array.newInstance(TestFieldProblemAbstract.class, 6); 234 allProblems[0] = new TestFieldProblem1<T>(field); 235 allProblems[1] = new TestFieldProblem2<T>(field); 236 allProblems[2] = new TestFieldProblem3<T>(field); 237 allProblems[3] = new TestFieldProblem4<T>(field); 238 allProblems[4] = new TestFieldProblem5<T>(field); 239 allProblems[5] = new TestFieldProblem6<T>(field); 240 for (TestFieldProblemAbstract<T> pb : allProblems) { 241 242 T previousValueError = null; 243 T previousTimeError = null; 244 for (int i = 4; i < 10; ++i) { 245 246 T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(FastMath.pow(2.0, -i)); 247 248 RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step); 249 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ); 250 integ.addStepHandler(handler); 251 FieldEventHandler<T>[] functions = pb.getEventsHandlers(); 252 for (int l = 0; l < functions.length; ++l) { 253 integ.addEventHandler(functions[l], 254 Double.POSITIVE_INFINITY, 1.0e-6 * step.getReal(), 1000); 255 } 256 Assert.assertEquals(functions.length, integ.getEventHandlers().size()); 257 FieldODEStateAndDerivative<T> stop = integ.integrate(new FieldExpandableODE<T>(pb), 258 pb.getInitialState(), 259 pb.getFinalTime()); 260 if (functions.length == 0) { 261 Assert.assertEquals(pb.getFinalTime().getReal(), stop.getTime().getReal(), epsilonT); 262 } 263 264 T error = handler.getMaximalValueError(); 265 if (i > 4) { 266 Assert.assertTrue(error.subtract(previousValueError.abs().multiply(safetyValueFactor)).getReal() < 0); 267 } 268 previousValueError = error; 269 270 T timeError = handler.getMaximalTimeError(); 271 if (i > 4) { 272 Assert.assertTrue(timeError.subtract(previousTimeError.abs().multiply(safetyTimeFactor)).getReal() <= 0); 273 } 274 previousTimeError = timeError; 275 276 integ.clearEventHandlers(); 277 Assert.assertEquals(0, integ.getEventHandlers().size()); 278 } 279 280 } 281 282 } 283 284 @Test 285 public abstract void testSmallStep(); 286 doTestSmallStep(Field<T> field, final double epsilonLast, final double epsilonMaxValue, final double epsilonMaxTime, final String name)287 protected <T extends RealFieldElement<T>> void doTestSmallStep(Field<T> field, 288 final double epsilonLast, 289 final double epsilonMaxValue, 290 final double epsilonMaxTime, 291 final String name) 292 throws DimensionMismatchException, NumberIsTooSmallException, 293 MaxCountExceededException, NoBracketingException { 294 295 TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field); 296 T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.001); 297 298 RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step); 299 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ); 300 integ.addStepHandler(handler); 301 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime()); 302 303 Assert.assertEquals(0, handler.getLastError().getReal(), epsilonLast); 304 Assert.assertEquals(0, handler.getMaximalValueError().getReal(), epsilonMaxValue); 305 Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), epsilonMaxTime); 306 Assert.assertEquals(name, integ.getName()); 307 308 } 309 310 @Test 311 public abstract void testBigStep(); 312 doTestBigStep(Field<T> field, final double belowLast, final double belowMaxValue, final double epsilonMaxTime, final String name)313 protected <T extends RealFieldElement<T>> void doTestBigStep(Field<T> field, 314 final double belowLast, 315 final double belowMaxValue, 316 final double epsilonMaxTime, 317 final String name) 318 throws DimensionMismatchException, NumberIsTooSmallException, 319 MaxCountExceededException, NoBracketingException { 320 321 TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field); 322 T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.2); 323 324 RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step); 325 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ); 326 integ.addStepHandler(handler); 327 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime()); 328 329 Assert.assertTrue(handler.getLastError().getReal() > belowLast); 330 Assert.assertTrue(handler.getMaximalValueError().getReal() > belowMaxValue); 331 Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), epsilonMaxTime); 332 Assert.assertEquals(name, integ.getName()); 333 334 } 335 336 @Test 337 public abstract void testBackward(); 338 doTestBackward(Field<T> field, final double epsilonLast, final double epsilonMaxValue, final double epsilonMaxTime, final String name)339 protected <T extends RealFieldElement<T>> void doTestBackward(Field<T> field, 340 final double epsilonLast, 341 final double epsilonMaxValue, 342 final double epsilonMaxTime, 343 final String name) 344 throws DimensionMismatchException, NumberIsTooSmallException, 345 MaxCountExceededException, NoBracketingException { 346 347 TestFieldProblem5<T> pb = new TestFieldProblem5<T>(field); 348 T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.001).abs(); 349 350 RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step); 351 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ); 352 integ.addStepHandler(handler); 353 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime()); 354 355 Assert.assertEquals(0, handler.getLastError().getReal(), epsilonLast); 356 Assert.assertEquals(0, handler.getMaximalValueError().getReal(), epsilonMaxValue); 357 Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), epsilonMaxTime); 358 Assert.assertEquals(name, integ.getName()); 359 360 } 361 362 @Test 363 public abstract void testKepler(); 364 doTestKepler(Field<T> field, double expectedMaxError, double epsilon)365 protected <T extends RealFieldElement<T>> void doTestKepler(Field<T> field, double expectedMaxError, double epsilon) 366 throws DimensionMismatchException, NumberIsTooSmallException, 367 MaxCountExceededException, NoBracketingException { 368 369 final TestFieldProblem3<T> pb = new TestFieldProblem3<T>(field, field.getZero().add(0.9)); 370 T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.0003); 371 372 RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step); 373 integ.addStepHandler(new KeplerHandler<T>(pb, expectedMaxError, epsilon)); 374 integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime()); 375 } 376 377 private static class KeplerHandler<T extends RealFieldElement<T>> implements FieldStepHandler<T> { 378 private T maxError; 379 private final TestFieldProblem3<T> pb; 380 private final double expectedMaxError; 381 private final double epsilon; KeplerHandler(TestFieldProblem3<T> pb, double expectedMaxError, double epsilon)382 public KeplerHandler(TestFieldProblem3<T> pb, double expectedMaxError, double epsilon) { 383 this.pb = pb; 384 this.expectedMaxError = expectedMaxError; 385 this.epsilon = epsilon; 386 maxError = pb.getField().getZero(); 387 } init(FieldODEStateAndDerivative<T> state0, T t)388 public void init(FieldODEStateAndDerivative<T> state0, T t) { 389 maxError = pb.getField().getZero(); 390 } handleStep(FieldStepInterpolator<T> interpolator, boolean isLast)391 public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast) 392 throws MaxCountExceededException { 393 394 FieldODEStateAndDerivative<T> current = interpolator.getCurrentState(); 395 T[] theoreticalY = pb.computeTheoreticalState(current.getTime()); 396 T dx = current.getState()[0].subtract(theoreticalY[0]); 397 T dy = current.getState()[1].subtract(theoreticalY[1]); 398 T error = dx.multiply(dx).add(dy.multiply(dy)); 399 if (error.subtract(maxError).getReal() > 0) { 400 maxError = error; 401 } 402 if (isLast) { 403 Assert.assertEquals(expectedMaxError, maxError.getReal(), epsilon); 404 } 405 } 406 } 407 408 @Test 409 public abstract void testStepSize(); 410 doTestStepSize(final Field<T> field, final double epsilon)411 protected <T extends RealFieldElement<T>> void doTestStepSize(final Field<T> field, final double epsilon) 412 throws DimensionMismatchException, NumberIsTooSmallException, 413 MaxCountExceededException, NoBracketingException { 414 final T step = field.getZero().add(1.23456); 415 RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step); 416 integ.addStepHandler(new FieldStepHandler<T>() { 417 public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast) { 418 if (! isLast) { 419 Assert.assertEquals(step.getReal(), 420 interpolator.getCurrentState().getTime().subtract(interpolator.getPreviousState().getTime()).getReal(), 421 epsilon); 422 } 423 } 424 public void init(FieldODEStateAndDerivative<T> s0, T t) { 425 } 426 }); 427 integ.integrate(new FieldExpandableODE<T>(new FirstOrderFieldDifferentialEquations<T>() { 428 public void init(T t0, T[] y0, T t) { 429 } 430 public T[] computeDerivatives(T t, T[] y) { 431 T[] dot = MathArrays.buildArray(t.getField(), 1); 432 dot[0] = t.getField().getOne(); 433 return dot; 434 } 435 public int getDimension() { 436 return 1; 437 } 438 }), new FieldODEState<T>(field.getZero(), MathArrays.buildArray(field, 1)), field.getZero().add(5.0)); 439 } 440 441 @Test 442 public abstract void testSingleStep(); 443 doTestSingleStep(final Field<T> field, final double epsilon)444 protected <T extends RealFieldElement<T>> void doTestSingleStep(final Field<T> field, final double epsilon) { 445 446 final TestFieldProblem3<T> pb = new TestFieldProblem3<T>(field, field.getZero().add(0.9)); 447 T h = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.0003); 448 449 RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, field.getZero().add(Double.NaN)); 450 T t = pb.getInitialState().getTime(); 451 T[] y = pb.getInitialState().getState(); 452 for (int i = 0; i < 100; ++i) { 453 y = integ.singleStep(pb, t, y, t.add(h)); 454 t = t.add(h); 455 } 456 T[] yth = pb.computeTheoreticalState(t); 457 T dx = y[0].subtract(yth[0]); 458 T dy = y[1].subtract(yth[1]); 459 T error = dx.multiply(dx).add(dy.multiply(dy)); 460 Assert.assertEquals(0.0, error.getReal(), epsilon); 461 } 462 463 @Test 464 public abstract void testTooLargeFirstStep(); 465 doTestTooLargeFirstStep(final Field<T> field)466 protected <T extends RealFieldElement<T>> void doTestTooLargeFirstStep(final Field<T> field) { 467 468 RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, field.getZero().add(0.5)); 469 final T t0 = field.getZero(); 470 final T[] y0 = MathArrays.buildArray(field, 1); 471 y0[0] = field.getOne(); 472 final T t = field.getZero().add(0.001); 473 FirstOrderFieldDifferentialEquations<T> equations = new FirstOrderFieldDifferentialEquations<T>() { 474 475 public int getDimension() { 476 return 1; 477 } 478 479 public void init(T t0, T[] y0, T t) { 480 } 481 482 public T[] computeDerivatives(T t, T[] y) { 483 Assert.assertTrue(t.getReal() >= FastMath.nextAfter(t0.getReal(), Double.NEGATIVE_INFINITY)); 484 Assert.assertTrue(t.getReal() <= FastMath.nextAfter(t.getReal(), Double.POSITIVE_INFINITY)); 485 T[] yDot = MathArrays.buildArray(field, 1); 486 yDot[0] = y[0].multiply(-100.0); 487 return yDot; 488 } 489 490 }; 491 492 integ.integrate(new FieldExpandableODE<T>(equations), new FieldODEState<T>(t0, y0), t); 493 494 } 495 496 @Test 497 public abstract void testUnstableDerivative(); 498 doTestUnstableDerivative(Field<T> field, double epsilon)499 protected <T extends RealFieldElement<T>> void doTestUnstableDerivative(Field<T> field, double epsilon) { 500 final StepFieldProblem<T> stepProblem = new StepFieldProblem<T>(field, 501 field.getZero().add(0.0), 502 field.getZero().add(1.0), 503 field.getZero().add(2.0)); 504 RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, field.getZero().add(0.3)); 505 integ.addEventHandler(stepProblem, 1.0, 1.0e-12, 1000); 506 FieldODEStateAndDerivative<T> result = integ.integrate(new FieldExpandableODE<T>(stepProblem), 507 new FieldODEState<T>(field.getZero(), MathArrays.buildArray(field, 1)), 508 field.getZero().add(10.0)); 509 Assert.assertEquals(8.0, result.getState()[0].getReal(), epsilon); 510 } 511 512 @Test 513 public abstract void testDerivativesConsistency(); 514 doTestDerivativesConsistency(final Field<T> field, double epsilon)515 protected <T extends RealFieldElement<T>> void doTestDerivativesConsistency(final Field<T> field, double epsilon) { 516 TestFieldProblem3<T> pb = new TestFieldProblem3<T>(field); 517 T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.001); 518 RungeKuttaFieldIntegrator<T> integ = createIntegrator(field, step); 519 StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 1.0e-10); 520 } 521 522 @Test 523 public abstract void testPartialDerivatives(); 524 doTestPartialDerivatives(final double epsilonY, final double[] epsilonPartials)525 protected void doTestPartialDerivatives(final double epsilonY, 526 final double[] epsilonPartials) { 527 528 // parameters indices 529 final int parameters = 5; 530 final int order = 1; 531 final int parOmega = 0; 532 final int parTO = 1; 533 final int parY00 = 2; 534 final int parY01 = 3; 535 final int parT = 4; 536 537 DerivativeStructure omega = new DerivativeStructure(parameters, order, parOmega, 1.3); 538 DerivativeStructure t0 = new DerivativeStructure(parameters, order, parTO, 1.3); 539 DerivativeStructure[] y0 = new DerivativeStructure[] { 540 new DerivativeStructure(parameters, order, parY00, 3.0), 541 new DerivativeStructure(parameters, order, parY01, 4.0) 542 }; 543 DerivativeStructure t = new DerivativeStructure(parameters, order, parT, 6.0); 544 SinCos sinCos = new SinCos(omega); 545 546 RungeKuttaFieldIntegrator<DerivativeStructure> integrator = 547 createIntegrator(omega.getField(), t.subtract(t0).multiply(0.001)); 548 FieldODEStateAndDerivative<DerivativeStructure> result = 549 integrator.integrate(new FieldExpandableODE<DerivativeStructure>(sinCos), 550 new FieldODEState<DerivativeStructure>(t0, y0), 551 t); 552 553 // check values 554 for (int i = 0; i < sinCos.getDimension(); ++i) { 555 Assert.assertEquals(sinCos.theoreticalY(t.getReal())[i], result.getState()[i].getValue(), epsilonY); 556 } 557 558 // check derivatives 559 final double[][] derivatives = sinCos.getDerivatives(t.getReal()); 560 for (int i = 0; i < sinCos.getDimension(); ++i) { 561 for (int parameter = 0; parameter < parameters; ++parameter) { 562 Assert.assertEquals(derivatives[i][parameter], 563 dYdP(result.getState()[i], parameter), 564 epsilonPartials[parameter]); 565 } 566 } 567 568 } 569 dYdP(final DerivativeStructure y, final int parameter)570 private double dYdP(final DerivativeStructure y, final int parameter) { 571 int[] orders = new int[y.getFreeParameters()]; 572 orders[parameter] = 1; 573 return y.getPartialDerivative(orders); 574 } 575 576 private static class SinCos implements FirstOrderFieldDifferentialEquations<DerivativeStructure> { 577 578 private final DerivativeStructure omega; 579 private DerivativeStructure r; 580 private DerivativeStructure alpha; 581 582 private double dRdY00; 583 private double dRdY01; 584 private double dAlphadOmega; 585 private double dAlphadT0; 586 private double dAlphadY00; 587 private double dAlphadY01; 588 SinCos(final DerivativeStructure omega)589 protected SinCos(final DerivativeStructure omega) { 590 this.omega = omega; 591 } 592 getDimension()593 public int getDimension() { 594 return 2; 595 } 596 init(final DerivativeStructure t0, final DerivativeStructure[] y0, final DerivativeStructure finalTime)597 public void init(final DerivativeStructure t0, final DerivativeStructure[] y0, 598 final DerivativeStructure finalTime) { 599 600 // theoretical solution is y(t) = { r * sin(omega * t + alpha), r * cos(omega * t + alpha) } 601 // so we retrieve alpha by identification from the initial state 602 final DerivativeStructure r2 = y0[0].multiply(y0[0]).add(y0[1].multiply(y0[1])); 603 604 this.r = r2.sqrt(); 605 this.dRdY00 = y0[0].divide(r).getReal(); 606 this.dRdY01 = y0[1].divide(r).getReal(); 607 608 this.alpha = y0[0].atan2(y0[1]).subtract(t0.multiply(omega)); 609 this.dAlphadOmega = -t0.getReal(); 610 this.dAlphadT0 = -omega.getReal(); 611 this.dAlphadY00 = y0[1].divide(r2).getReal(); 612 this.dAlphadY01 = y0[0].negate().divide(r2).getReal(); 613 614 } 615 computeDerivatives(final DerivativeStructure t, final DerivativeStructure[] y)616 public DerivativeStructure[] computeDerivatives(final DerivativeStructure t, final DerivativeStructure[] y) { 617 return new DerivativeStructure[] { 618 omega.multiply(y[1]), 619 omega.multiply(y[0]).negate() 620 }; 621 } 622 theoreticalY(final double t)623 public double[] theoreticalY(final double t) { 624 final double theta = omega.getReal() * t + alpha.getReal(); 625 return new double[] { 626 r.getReal() * FastMath.sin(theta), r.getReal() * FastMath.cos(theta) 627 }; 628 } 629 getDerivatives(final double t)630 public double[][] getDerivatives(final double t) { 631 632 // intermediate angle and state 633 final double theta = omega.getReal() * t + alpha.getReal(); 634 final double sin = FastMath.sin(theta); 635 final double cos = FastMath.cos(theta); 636 final double y0 = r.getReal() * sin; 637 final double y1 = r.getReal() * cos; 638 639 // partial derivatives of the state first component 640 final double dY0dOmega = y1 * (t + dAlphadOmega); 641 final double dY0dT0 = y1 * dAlphadT0; 642 final double dY0dY00 = dRdY00 * sin + y1 * dAlphadY00; 643 final double dY0dY01 = dRdY01 * sin + y1 * dAlphadY01; 644 final double dY0dT = y1 * omega.getReal(); 645 646 // partial derivatives of the state second component 647 final double dY1dOmega = - y0 * (t + dAlphadOmega); 648 final double dY1dT0 = - y0 * dAlphadT0; 649 final double dY1dY00 = dRdY00 * cos - y0 * dAlphadY00; 650 final double dY1dY01 = dRdY01 * cos - y0 * dAlphadY01; 651 final double dY1dT = - y0 * omega.getReal(); 652 653 return new double[][] { 654 { dY0dOmega, dY0dT0, dY0dY00, dY0dY01, dY0dT }, 655 { dY1dOmega, dY1dT0, dY1dY00, dY1dY01, dY1dT } 656 }; 657 658 } 659 660 } 661 662 } 663