View Javadoc
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.math4.legacy.ode;
19  
20  import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
21  import org.apache.commons.math4.legacy.ode.sampling.StepHandler;
22  import org.apache.commons.math4.legacy.ode.sampling.StepInterpolator;
23  import org.apache.commons.math4.core.jdkmath.JdkMath;
24  
25  /**
26   * This class is used to handle steps for the test problems
27   * integrated during the junit tests for the ODE integrators.
28   */
29  public class TestProblemHandler
30    implements StepHandler {
31  
32    /** Associated problem. */
33    private TestProblemAbstract problem;
34  
35    /** Maximal errors encountered during the integration. */
36    private double maxValueError;
37    private double maxTimeError;
38  
39    /** Error at the end of the integration. */
40    private double lastError;
41  
42    /** Time at the end of integration. */
43    private double lastTime;
44  
45    /** ODE solver used. */
46    private ODEIntegrator integrator;
47  
48    /** Expected start for step. */
49    private double expectedStepStart;
50  
51    /**
52     * Simple constructor.
53     * @param problem problem for which steps should be handled
54     * @param integrator ODE solver used
55     */
56    public TestProblemHandler(TestProblemAbstract problem, ODEIntegrator integrator) {
57      this.problem = problem;
58      this.integrator = integrator;
59      maxValueError = 0;
60      maxTimeError  = 0;
61      lastError     = 0;
62      expectedStepStart = Double.NaN;
63    }
64  
65    @Override
66  public void init(double t0, double[] y0, double t) {
67      maxValueError = 0;
68      maxTimeError  = 0;
69      lastError     = 0;
70      expectedStepStart = Double.NaN;
71    }
72  
73    @Override
74  public void handleStep(StepInterpolator interpolator, boolean isLast) throws MaxCountExceededException {
75  
76      double start = integrator.getCurrentStepStart();
77      if (JdkMath.abs((start - problem.getInitialTime()) / integrator.getCurrentSignedStepsize()) > 0.001) {
78          // multistep integrators do not handle the first steps themselves
79          // so we have to make sure the integrator we look at has really started its work
80          if (!Double.isNaN(expectedStepStart)) {
81              // the step should either start at the end of the integrator step
82              // or at an event if the step is split into several substeps
83              double stepError = JdkMath.max(maxTimeError, JdkMath.abs(start - expectedStepStart));
84              for (double eventTime : problem.getTheoreticalEventsTimes()) {
85                  stepError = JdkMath.min(stepError, JdkMath.abs(start - eventTime));
86              }
87              maxTimeError = JdkMath.max(maxTimeError, stepError);
88          }
89          expectedStepStart = start + integrator.getCurrentSignedStepsize();
90      }
91  
92  
93      double pT = interpolator.getPreviousTime();
94      double cT = interpolator.getCurrentTime();
95      double[] errorScale = problem.getErrorScale();
96  
97      // store the error at the last step
98      if (isLast) {
99          double[] interpolatedY = interpolator.getInterpolatedState();
100         double[] theoreticalY  = problem.computeTheoreticalState(cT);
101         for (int i = 0; i < interpolatedY.length; ++i) {
102             double error = JdkMath.abs(interpolatedY[i] - theoreticalY[i]);
103             lastError = JdkMath.max(error, lastError);
104         }
105         lastTime = cT;
106     }
107     // walk through the step
108     for (int k = 0; k <= 20; ++k) {
109 
110         double time = pT + (k * (cT - pT)) / 20;
111         interpolator.setInterpolatedTime(time);
112         double[] interpolatedY = interpolator.getInterpolatedState();
113         double[] theoreticalY  = problem.computeTheoreticalState(interpolator.getInterpolatedTime());
114 
115         // update the errors
116         for (int i = 0; i < interpolatedY.length; ++i) {
117             double error = errorScale[i] * JdkMath.abs(interpolatedY[i] - theoreticalY[i]);
118             maxValueError = JdkMath.max(error, maxValueError);
119         }
120     }
121   }
122 
123   /**
124    * Get the maximal value error encountered during integration.
125    * @return maximal value error
126    */
127   public double getMaximalValueError() {
128     return maxValueError;
129   }
130 
131   /**
132    * Get the maximal time error encountered during integration.
133    * @return maximal time error
134    */
135   public double getMaximalTimeError() {
136     return maxTimeError;
137   }
138 
139 
140   public int getCalls() {
141       return problem.getCalls();
142   }
143 
144   /**
145    * Get the error at the end of the integration.
146    * @return error at the end of the integration
147    */
148   public double getLastError() {
149     return lastError;
150   }
151 
152   /**
153    * Get the time at the end of the integration.
154    * @return time at the end of the integration.
155    */
156   public double getLastTime() {
157     return lastTime;
158   }
159 }