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