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.core.RealFieldElement;
21  import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
22  import org.apache.commons.math4.legacy.ode.sampling.FieldStepHandler;
23  import org.apache.commons.math4.legacy.ode.sampling.FieldStepInterpolator;
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   * @param <T> the type of the field elements
29   */
30  public class TestFieldProblemHandler<T extends RealFieldElement<T>>
31      implements FieldStepHandler<T> {
32  
33      /** Associated problem. */
34      private TestFieldProblemAbstract<T> problem;
35  
36      /** Maximal errors encountered during the integration. */
37      private T maxValueError;
38      private T maxTimeError;
39  
40      /** Error at the end of the integration. */
41      private T lastError;
42  
43      /** Time at the end of integration. */
44      private T lastTime;
45  
46      /** ODE solver used. */
47      private FirstOrderFieldIntegrator<T> integrator;
48  
49      /** Expected start for step. */
50      private T expectedStepStart;
51  
52      /**
53       * Simple constructor.
54       * @param problem problem for which steps should be handled
55       * @param integrator ODE solver used
56       */
57      public TestFieldProblemHandler(TestFieldProblemAbstract<T> problem, FirstOrderFieldIntegrator<T> integrator) {
58          this.problem      = problem;
59          this.integrator   = integrator;
60          maxValueError     = problem.getField().getZero();
61          maxTimeError      = problem.getField().getZero();
62          lastError         = problem.getField().getZero();
63          expectedStepStart = null;
64      }
65  
66      @Override
67      public void init(FieldODEStateAndDerivative<T> state0, T t) {
68          maxValueError     = problem.getField().getZero();
69          maxTimeError      = problem.getField().getZero();
70          lastError         = problem.getField().getZero();
71          expectedStepStart = null;
72      }
73  
74      @Override
75      public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast) throws MaxCountExceededException {
76  
77          T start = integrator.getCurrentStepStart().getTime();
78          if (start.subtract(problem.getInitialState().getTime()).divide(integrator.getCurrentSignedStepsize()).abs().getReal() > 0.001) {
79              // multistep integrators do not handle the first steps themselves
80              // so we have to make sure the integrator we look at has really started its work
81              if (expectedStepStart != null) {
82                  // the step should either start at the end of the integrator step
83                  // or at an event if the step is split into several substeps
84                  T stepError = RealFieldElement.max(maxTimeError, start.subtract(expectedStepStart).abs());
85                  for (T eventTime : problem.getTheoreticalEventsTimes()) {
86                      stepError = RealFieldElement.min(stepError, start.subtract(eventTime).abs());
87                  }
88                  maxTimeError = RealFieldElement.max(maxTimeError, stepError);
89              }
90              expectedStepStart = start.add(integrator.getCurrentSignedStepsize());
91          }
92  
93          T pT = interpolator.getPreviousState().getTime();
94          T cT = interpolator.getCurrentState().getTime();
95          T[] errorScale = problem.getErrorScale();
96  
97          // store the error at the last step
98          if (isLast) {
99              T[] interpolatedY = interpolator.getInterpolatedState(cT).getState();
100             T[] theoreticalY  = problem.computeTheoreticalState(cT);
101             for (int i = 0; i < interpolatedY.length; ++i) {
102                 T error = interpolatedY[i].subtract(theoreticalY[i]).abs();
103                 lastError = RealFieldElement.max(error, lastError);
104             }
105             lastTime = cT;
106         }
107 
108         // walk through the step
109         for (int k = 0; k <= 20; ++k) {
110 
111             T time = pT.add(cT.subtract(pT).multiply(k).divide(20));
112             T[] interpolatedY = interpolator.getInterpolatedState(time).getState();
113             T[] theoreticalY  = problem.computeTheoreticalState(time);
114 
115             // update the errors
116             for (int i = 0; i < interpolatedY.length; ++i) {
117                 T error = errorScale[i].multiply(interpolatedY[i].subtract(theoreticalY[i]).abs());
118                 maxValueError = RealFieldElement.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 T getMaximalValueError() {
128         return maxValueError;
129     }
130 
131     /**
132      * Get the maximal time error encountered during integration.
133      * @return maximal time error
134      */
135     public T getMaximalTimeError() {
136         return maxTimeError;
137     }
138 
139     /**
140      * Get the error at the end of the integration.
141      * @return error at the end of the integration
142      */
143     public T getLastError() {
144         return lastError;
145     }
146 
147     /**
148      * Get the time at the end of the integration.
149      * @return time at the end of the integration.
150      */
151     public T getLastTime() {
152         return lastTime;
153     }
154 }