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  package org.apache.commons.math.ode.sampling;
18  
19  import static org.junit.Assert.assertEquals;
20  
21  import org.apache.commons.math.ode.DerivativeException;
22  import org.apache.commons.math.ode.FirstOrderIntegrator;
23  import org.apache.commons.math.ode.IntegratorException;
24  import org.apache.commons.math.ode.TestProblemAbstract;
25  
26  public class StepInterpolatorTestUtils {
27  
28      public static void checkDerivativesConsistency(final FirstOrderIntegrator integrator,
29                                                     final TestProblemAbstract problem,
30                                                     final double threshold)
31          throws DerivativeException, IntegratorException {
32          integrator.addStepHandler(new StepHandler() {
33  
34              public boolean requiresDenseOutput() {
35                  return true;
36              }
37  
38              public void handleStep(StepInterpolator interpolator, boolean isLast)
39                  throws DerivativeException {
40  
41                  final double h = 0.001 * (interpolator.getCurrentTime() - interpolator.getPreviousTime());
42                  final double t = interpolator.getCurrentTime() - 300 * h;
43  
44                  if (Math.abs(h) < 10 * Math.ulp(t)) {
45                      return;
46                  }
47  
48                  interpolator.setInterpolatedTime(t - 4 * h);
49                  final double[] yM4h = interpolator.getInterpolatedState().clone();
50                  interpolator.setInterpolatedTime(t - 3 * h);
51                  final double[] yM3h = interpolator.getInterpolatedState().clone();
52                  interpolator.setInterpolatedTime(t - 2 * h);
53                  final double[] yM2h = interpolator.getInterpolatedState().clone();
54                  interpolator.setInterpolatedTime(t - h);
55                  final double[] yM1h = interpolator.getInterpolatedState().clone();
56                  interpolator.setInterpolatedTime(t + h);
57                  final double[] yP1h = interpolator.getInterpolatedState().clone();
58                  interpolator.setInterpolatedTime(t + 2 * h);
59                  final double[] yP2h = interpolator.getInterpolatedState().clone();
60                  interpolator.setInterpolatedTime(t + 3 * h);
61                  final double[] yP3h = interpolator.getInterpolatedState().clone();
62                  interpolator.setInterpolatedTime(t + 4 * h);
63                  final double[] yP4h = interpolator.getInterpolatedState().clone();
64  
65                  interpolator.setInterpolatedTime(t);
66                  final double[] yDot = interpolator.getInterpolatedDerivatives();
67  
68                  for (int i = 0; i < yDot.length; ++i) {
69                      final double approYDot = ( -3 * (yP4h[i] - yM4h[i]) +
70                                                 32 * (yP3h[i] - yM3h[i]) +
71                                               -168 * (yP2h[i] - yM2h[i]) +
72                                                672 * (yP1h[i] - yM1h[i])) / (840 * h);
73                      assertEquals(approYDot, yDot[i], threshold);
74                  }
75  
76              }
77  
78              public void reset() {
79              }
80  
81          });
82  
83          integrator.integrate(problem,
84                               problem.getInitialTime(), problem.getInitialState(),
85                               problem.getFinalTime(), new double[problem.getDimension()]);
86  
87      }
88  }
89