1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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
27
28
29
30 public class TestFieldProblemHandler<T extends RealFieldElement<T>>
31 implements FieldStepHandler<T> {
32
33
34 private TestFieldProblemAbstract<T> problem;
35
36
37 private T maxValueError;
38 private T maxTimeError;
39
40
41 private T lastError;
42
43
44 private T lastTime;
45
46
47 private FirstOrderFieldIntegrator<T> integrator;
48
49
50 private T expectedStepStart;
51
52
53
54
55
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
80
81 if (expectedStepStart != null) {
82
83
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
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
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
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
125
126
127 public T getMaximalValueError() {
128 return maxValueError;
129 }
130
131
132
133
134
135 public T getMaximalTimeError() {
136 return maxTimeError;
137 }
138
139
140
141
142
143 public T getLastError() {
144 return lastError;
145 }
146
147
148
149
150
151 public T getLastTime() {
152 return lastTime;
153 }
154 }