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.nonstiff;
19
20
21 import org.apache.commons.math4.legacy.core.Field;
22 import org.apache.commons.math4.legacy.core.RealFieldElement;
23 import org.apache.commons.math4.legacy.exception.MathIllegalStateException;
24 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
25 import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
26 import org.apache.commons.math4.legacy.ode.AbstractFieldIntegrator;
27 import org.apache.commons.math4.legacy.ode.FieldExpandableODE;
28 import org.apache.commons.math4.legacy.ode.FieldODEState;
29 import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative;
30 import org.apache.commons.math4.legacy.ode.FirstOrderFieldIntegrator;
31 import org.apache.commons.math4.legacy.ode.MultistepFieldIntegrator;
32 import org.apache.commons.math4.legacy.ode.TestFieldProblem1;
33 import org.apache.commons.math4.legacy.ode.TestFieldProblem5;
34 import org.apache.commons.math4.legacy.ode.TestFieldProblem6;
35 import org.apache.commons.math4.legacy.ode.TestFieldProblemAbstract;
36 import org.apache.commons.math4.legacy.ode.TestFieldProblemHandler;
37 import org.apache.commons.math4.legacy.ode.sampling.FieldStepHandler;
38 import org.apache.commons.math4.legacy.ode.sampling.FieldStepInterpolator;
39 import org.apache.commons.math4.core.jdkmath.JdkMath;
40 import org.junit.Assert;
41 import org.junit.Test;
42
43 public abstract class AdamsFieldIntegratorAbstractTest {
44
45 protected abstract <T extends RealFieldElement<T>> AdamsFieldIntegrator<T>
46 createIntegrator(Field<T> field, int nSteps, double minStep, double maxStep,
47 double scalAbsoluteTolerance, double scalRelativeTolerance);
48
49 protected abstract <T extends RealFieldElement<T>> AdamsFieldIntegrator<T>
50 createIntegrator(Field<T> field, int nSteps, double minStep, double maxStep,
51 double[] vecAbsoluteTolerance, double[] vecRelativeTolerance);
52
53 @Test(expected=NumberIsTooSmallException.class)
54 public abstract void testMinStep();
55
56 protected <T extends RealFieldElement<T>> void doDimensionCheck(final Field<T> field) {
57 TestFieldProblem1<T> pb = new TestFieldProblem1<>(field);
58
59 double minStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.1).getReal();
60 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
61 double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
62 double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
63
64 FirstOrderFieldIntegrator<T> integ = createIntegrator(field, 4, minStep, maxStep,
65 vecAbsoluteTolerance,
66 vecRelativeTolerance);
67 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<>(pb, integ);
68 integ.addStepHandler(handler);
69 integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
70 }
71
72 @Test
73 public abstract void testIncreasingTolerance();
74
75 protected <T extends RealFieldElement<T>> void doTestIncreasingTolerance(final Field<T> field,
76 double ratioMin, double ratioMax) {
77
78 int previousCalls = Integer.MAX_VALUE;
79 for (int i = -12; i < -2; ++i) {
80 TestFieldProblem1<T> pb = new TestFieldProblem1<>(field);
81 double minStep = 0;
82 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
83 double scalAbsoluteTolerance = JdkMath.pow(10.0, i);
84 double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
85
86 FirstOrderFieldIntegrator<T> integ = createIntegrator(field, 4, minStep, maxStep,
87 scalAbsoluteTolerance,
88 scalRelativeTolerance);
89 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<>(pb, integ);
90 integ.addStepHandler(handler);
91 integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
92
93 Assert.assertTrue(handler.getMaximalValueError().getReal() > ratioMin * scalAbsoluteTolerance);
94 Assert.assertTrue(handler.getMaximalValueError().getReal() < ratioMax * scalAbsoluteTolerance);
95
96 int calls = pb.getCalls();
97 Assert.assertEquals(integ.getEvaluations(), calls);
98 Assert.assertTrue(calls <= previousCalls);
99 previousCalls = calls;
100 }
101 }
102
103 @Test(expected = MaxCountExceededException.class)
104 public abstract void exceedMaxEvaluations();
105
106 protected <T extends RealFieldElement<T>> void doExceedMaxEvaluations(final Field<T> field, final int max) {
107
108 TestFieldProblem1<T> pb = new TestFieldProblem1<>(field);
109 double range = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
110
111 FirstOrderFieldIntegrator<T> integ = createIntegrator(field, 2, 0, range, 1.0e-12, 1.0e-12);
112 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<>(pb, integ);
113 integ.addStepHandler(handler);
114 integ.setMaxEvaluations(max);
115 integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
116 }
117
118 @Test
119 public abstract void backward();
120
121 protected <T extends RealFieldElement<T>> void doBackward(final Field<T> field,
122 final double epsilonLast,
123 final double epsilonMaxValue,
124 final double epsilonMaxTime,
125 final String name) {
126
127 TestFieldProblem5<T> pb = new TestFieldProblem5<>(field);
128 double range = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
129
130 AdamsFieldIntegrator<T> integ = createIntegrator(field, 4, 0, range, 1.0e-12, 1.0e-12);
131 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<>(pb, integ);
132 integ.addStepHandler(handler);
133 integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
134
135 Assert.assertEquals(0.0, handler.getLastError().getReal(), epsilonLast);
136 Assert.assertEquals(0.0, handler.getMaximalValueError().getReal(), epsilonMaxValue);
137 Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), epsilonMaxTime);
138 Assert.assertEquals(name, integ.getName());
139 }
140
141 @Test
142 public abstract void polynomial();
143
144 protected <T extends RealFieldElement<T>> void doPolynomial(final Field<T> field,
145 final int nLimit,
146 final double epsilonBad,
147 final double epsilonGood) {
148 TestFieldProblem6<T> pb = new TestFieldProblem6<>(field);
149 double range = pb.getFinalTime().subtract(pb.getInitialState().getTime()).abs().getReal();
150
151 for (int nSteps = 2; nSteps < 8; ++nSteps) {
152 AdamsFieldIntegrator<T> integ = createIntegrator(field, nSteps, 1.0e-6 * range, 0.1 * range, 1.0e-4, 1.0e-4);
153 integ.setStarterIntegrator(new PerfectStarter<>(pb, nSteps));
154 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<>(pb, integ);
155 integ.addStepHandler(handler);
156 integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
157 if (nSteps < nLimit) {
158 Assert.assertTrue(handler.getMaximalValueError().getReal() > epsilonBad);
159 } else {
160 Assert.assertTrue(handler.getMaximalValueError().getReal() < epsilonGood);
161 }
162 }
163 }
164
165 @Test(expected=MathIllegalStateException.class)
166 public abstract void testStartFailure();
167
168 protected <T extends RealFieldElement<T>> void doTestStartFailure(final Field<T> field) {
169 TestFieldProblem1<T> pb = new TestFieldProblem1<>(field);
170 double minStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.0001).getReal();
171 double maxStep = pb.getFinalTime().subtract(pb.getInitialState().getTime()).getReal();
172 double scalAbsoluteTolerance = 1.0e-6;
173 double scalRelativeTolerance = 1.0e-7;
174
175 MultistepFieldIntegrator<T> integ = createIntegrator(field, 6, minStep, maxStep,
176 scalAbsoluteTolerance,
177 scalRelativeTolerance);
178 integ.setStarterIntegrator(new DormandPrince853FieldIntegrator<>(field, maxStep * 0.5, maxStep, 0.1, 0.1));
179 TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<>(pb, integ);
180 integ.addStepHandler(handler);
181 integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
182 }
183
184 private static final class PerfectStarter<T extends RealFieldElement<T>> extends AbstractFieldIntegrator<T> {
185
186 private final PerfectInterpolator<T> interpolator;
187 private final int nbSteps;
188
189 PerfectStarter(final TestFieldProblemAbstract<T> problem, final int nbSteps) {
190 super(problem.getField(), "perfect-starter");
191 this.interpolator = new PerfectInterpolator<>(problem);
192 this.nbSteps = nbSteps;
193 }
194
195 @Override
196 public FieldODEStateAndDerivative<T> integrate(FieldExpandableODE<T> equations,
197 FieldODEState<T> initialState, T finalTime) {
198 T tStart = initialState.getTime().add(finalTime.subtract(initialState.getTime()).multiply(0.01));
199 getEvaluationsCounter().increment(nbSteps);
200 interpolator.setCurrentTime(initialState.getTime());
201 for (int i = 0; i < nbSteps; ++i) {
202 T tK = initialState.getTime().multiply(nbSteps - 1 - (i + 1)).add(tStart.multiply(i + 1)).divide(nbSteps - 1);
203 interpolator.setPreviousTime(interpolator.getCurrentTime());
204 interpolator.setCurrentTime(tK);
205 for (FieldStepHandler<T> handler : getStepHandlers()) {
206 handler.handleStep(interpolator, i == nbSteps - 1);
207 }
208 }
209 return interpolator.getInterpolatedState(tStart);
210 }
211 }
212
213 private static final class PerfectInterpolator<T extends RealFieldElement<T>> implements FieldStepInterpolator<T> {
214 private final TestFieldProblemAbstract<T> problem;
215 private T previousTime;
216 private T currentTime;
217
218 PerfectInterpolator(final TestFieldProblemAbstract<T> problem) {
219 this.problem = problem;
220 }
221
222 public void setPreviousTime(T previousTime) {
223 this.previousTime = previousTime;
224 }
225
226 public void setCurrentTime(T currentTime) {
227 this.currentTime = currentTime;
228 }
229
230 public T getCurrentTime() {
231 return currentTime;
232 }
233
234 @Override
235 public boolean isForward() {
236 return problem.getFinalTime().subtract(problem.getInitialState().getTime()).getReal() >= 0;
237 }
238
239 @Override
240 public FieldODEStateAndDerivative<T> getPreviousState() {
241 return getInterpolatedState(previousTime);
242 }
243
244 @Override
245 public FieldODEStateAndDerivative<T> getCurrentState() {
246 return getInterpolatedState(currentTime);
247 }
248
249 @Override
250 public FieldODEStateAndDerivative<T> getInterpolatedState(T time) {
251 T[] y = problem.computeTheoreticalState(time);
252 T[] yDot = problem.computeDerivatives(time, y);
253 return new FieldODEStateAndDerivative<>(time, y, yDot);
254 }
255 }
256 }