001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    package org.apache.commons.math.ode.sampling;
018    
019    import static org.junit.Assert.assertEquals;
020    
021    import org.apache.commons.math.ode.DerivativeException;
022    import org.apache.commons.math.ode.FirstOrderIntegrator;
023    import org.apache.commons.math.ode.IntegratorException;
024    import org.apache.commons.math.ode.TestProblemAbstract;
025    
026    public class StepInterpolatorTestUtils {
027    
028        public static void checkDerivativesConsistency(final FirstOrderIntegrator integrator,
029                                                       final TestProblemAbstract problem,
030                                                       final double threshold)
031            throws DerivativeException, IntegratorException {
032            integrator.addStepHandler(new StepHandler() {
033    
034                public boolean requiresDenseOutput() {
035                    return true;
036                }
037    
038                public void handleStep(StepInterpolator interpolator, boolean isLast)
039                    throws DerivativeException {
040    
041                    final double h = 0.001 * (interpolator.getCurrentTime() - interpolator.getPreviousTime());
042                    final double t = interpolator.getCurrentTime() - 300 * h;
043    
044                    if (Math.abs(h) < 10 * Math.ulp(t)) {
045                        return;
046                    }
047    
048                    interpolator.setInterpolatedTime(t - 4 * h);
049                    final double[] yM4h = interpolator.getInterpolatedState().clone();
050                    interpolator.setInterpolatedTime(t - 3 * h);
051                    final double[] yM3h = interpolator.getInterpolatedState().clone();
052                    interpolator.setInterpolatedTime(t - 2 * h);
053                    final double[] yM2h = interpolator.getInterpolatedState().clone();
054                    interpolator.setInterpolatedTime(t - h);
055                    final double[] yM1h = interpolator.getInterpolatedState().clone();
056                    interpolator.setInterpolatedTime(t + h);
057                    final double[] yP1h = interpolator.getInterpolatedState().clone();
058                    interpolator.setInterpolatedTime(t + 2 * h);
059                    final double[] yP2h = interpolator.getInterpolatedState().clone();
060                    interpolator.setInterpolatedTime(t + 3 * h);
061                    final double[] yP3h = interpolator.getInterpolatedState().clone();
062                    interpolator.setInterpolatedTime(t + 4 * h);
063                    final double[] yP4h = interpolator.getInterpolatedState().clone();
064    
065                    interpolator.setInterpolatedTime(t);
066                    final double[] yDot = interpolator.getInterpolatedDerivatives();
067    
068                    for (int i = 0; i < yDot.length; ++i) {
069                        final double approYDot = ( -3 * (yP4h[i] - yM4h[i]) +
070                                                   32 * (yP3h[i] - yM3h[i]) +
071                                                 -168 * (yP2h[i] - yM2h[i]) +
072                                                  672 * (yP1h[i] - yM1h[i])) / (840 * h);
073                        assertEquals(approYDot, yDot[i], threshold);
074                    }
075    
076                }
077    
078                public void reset() {
079                }
080    
081            });
082    
083            integrator.integrate(problem,
084                                 problem.getInitialTime(), problem.getInitialState(),
085                                 problem.getFinalTime(), new double[problem.getDimension()]);
086    
087        }
088    }
089