source: src/main/java/agents/anac/y2019/harddealer/math3/ode/MultistepFieldIntegrator.java

Last change on this file was 204, checked in by Katsuhide Fujita, 5 years ago

Fixed errors of ANAC2019 agents

  • Property svn:executable set to *
File size: 18.5 KB
Line 
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
18package agents.anac.y2019.harddealer.math3.ode;
19
20import agents.anac.y2019.harddealer.math3.Field;
21import agents.anac.y2019.harddealer.math3.RealFieldElement;
22import agents.anac.y2019.harddealer.math3.exception.DimensionMismatchException;
23import agents.anac.y2019.harddealer.math3.exception.MathIllegalStateException;
24import agents.anac.y2019.harddealer.math3.exception.MaxCountExceededException;
25import agents.anac.y2019.harddealer.math3.exception.NoBracketingException;
26import agents.anac.y2019.harddealer.math3.exception.NumberIsTooSmallException;
27import agents.anac.y2019.harddealer.math3.exception.util.LocalizedFormats;
28import agents.anac.y2019.harddealer.math3.linear.Array2DRowFieldMatrix;
29import agents.anac.y2019.harddealer.math3.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
30import agents.anac.y2019.harddealer.math3.ode.nonstiff.DormandPrince853FieldIntegrator;
31import agents.anac.y2019.harddealer.math3.ode.sampling.FieldStepHandler;
32import agents.anac.y2019.harddealer.math3.ode.sampling.FieldStepInterpolator;
33import agents.anac.y2019.harddealer.math3.util.FastMath;
34import agents.anac.y2019.harddealer.math3.util.MathArrays;
35import agents.anac.y2019.harddealer.math3.util.MathUtils;
36
37/**
38 * This class is the base class for multistep integrators for Ordinary
39 * Differential Equations.
40 * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
41 * <pre>
42 * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
43 * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
44 * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
45 * ...
46 * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
47 * </pre></p>
48 * <p>Rather than storing several previous steps separately, this implementation uses
49 * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
50 * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
51 * <pre>
52 * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
53 * </pre>
54 * (we omit the k index in the notation for clarity)</p>
55 * <p>
56 * Multistep integrators with Nordsieck representation are highly sensitive to
57 * large step changes because when the step is multiplied by factor a, the
58 * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
59 * and the last components are the least accurate ones. The default max growth
60 * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
61 * </p>
62 *
63 * @see agents.anac.y2019.harddealer.math3.ode.nonstiff.AdamsBashforthFieldIntegrator
64 * @see agents.anac.y2019.harddealer.math3.ode.nonstiff.AdamsMoultonFieldIntegrator
65 * @param <T> the type of the field elements
66 * @since 3.6
67 */
68public abstract class MultistepFieldIntegrator<T extends RealFieldElement<T>>
69 extends AdaptiveStepsizeFieldIntegrator<T> {
70
71 /** First scaled derivative (h y'). */
72 protected T[] scaled;
73
74 /** Nordsieck matrix of the higher scaled derivatives.
75 * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y<sup>(k)</sup>)</p>
76 */
77 protected Array2DRowFieldMatrix<T> nordsieck;
78
79 /** Starter integrator. */
80 private FirstOrderFieldIntegrator<T> starter;
81
82 /** Number of steps of the multistep method (excluding the one being computed). */
83 private final int nSteps;
84
85 /** Stepsize control exponent. */
86 private double exp;
87
88 /** Safety factor for stepsize control. */
89 private double safety;
90
91 /** Minimal reduction factor for stepsize control. */
92 private double minReduction;
93
94 /** Maximal growth factor for stepsize control. */
95 private double maxGrowth;
96
97 /**
98 * Build a multistep integrator with the given stepsize bounds.
99 * <p>The default starter integrator is set to the {@link
100 * DormandPrince853FieldIntegrator Dormand-Prince 8(5,3)} integrator with
101 * some defaults settings.</p>
102 * <p>
103 * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
104 * </p>
105 * @param field field to which the time and state vector elements belong
106 * @param name name of the method
107 * @param nSteps number of steps of the multistep method
108 * (excluding the one being computed)
109 * @param order order of the method
110 * @param minStep minimal step (must be positive even for backward
111 * integration), the last step can be smaller than this
112 * @param maxStep maximal step (must be positive even for backward
113 * integration)
114 * @param scalAbsoluteTolerance allowed absolute error
115 * @param scalRelativeTolerance allowed relative error
116 * @exception NumberIsTooSmallException if number of steps is smaller than 2
117 */
118 protected MultistepFieldIntegrator(final Field<T> field, final String name,
119 final int nSteps, final int order,
120 final double minStep, final double maxStep,
121 final double scalAbsoluteTolerance,
122 final double scalRelativeTolerance)
123 throws NumberIsTooSmallException {
124
125 super(field, name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
126
127 if (nSteps < 2) {
128 throw new NumberIsTooSmallException(
129 LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
130 nSteps, 2, true);
131 }
132
133 starter = new DormandPrince853FieldIntegrator<T>(field, minStep, maxStep,
134 scalAbsoluteTolerance,
135 scalRelativeTolerance);
136 this.nSteps = nSteps;
137
138 exp = -1.0 / order;
139
140 // set the default values of the algorithm control parameters
141 setSafety(0.9);
142 setMinReduction(0.2);
143 setMaxGrowth(FastMath.pow(2.0, -exp));
144
145 }
146
147 /**
148 * Build a multistep integrator with the given stepsize bounds.
149 * <p>The default starter integrator is set to the {@link
150 * DormandPrince853FieldIntegrator Dormand-Prince 8(5,3)} integrator with
151 * some defaults settings.</p>
152 * <p>
153 * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
154 * </p>
155 * @param field field to which the time and state vector elements belong
156 * @param name name of the method
157 * @param nSteps number of steps of the multistep method
158 * (excluding the one being computed)
159 * @param order order of the method
160 * @param minStep minimal step (must be positive even for backward
161 * integration), the last step can be smaller than this
162 * @param maxStep maximal step (must be positive even for backward
163 * integration)
164 * @param vecAbsoluteTolerance allowed absolute error
165 * @param vecRelativeTolerance allowed relative error
166 */
167 protected MultistepFieldIntegrator(final Field<T> field, final String name, final int nSteps,
168 final int order,
169 final double minStep, final double maxStep,
170 final double[] vecAbsoluteTolerance,
171 final double[] vecRelativeTolerance) {
172 super(field, name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
173 starter = new DormandPrince853FieldIntegrator<T>(field, minStep, maxStep,
174 vecAbsoluteTolerance,
175 vecRelativeTolerance);
176 this.nSteps = nSteps;
177
178 exp = -1.0 / order;
179
180 // set the default values of the algorithm control parameters
181 setSafety(0.9);
182 setMinReduction(0.2);
183 setMaxGrowth(FastMath.pow(2.0, -exp));
184
185 }
186
187 /**
188 * Get the starter integrator.
189 * @return starter integrator
190 */
191 public FirstOrderFieldIntegrator<T> getStarterIntegrator() {
192 return starter;
193 }
194
195 /**
196 * Set the starter integrator.
197 * <p>The various step and event handlers for this starter integrator
198 * will be managed automatically by the multi-step integrator. Any
199 * user configuration for these elements will be cleared before use.</p>
200 * @param starterIntegrator starter integrator
201 */
202 public void setStarterIntegrator(FirstOrderFieldIntegrator<T> starterIntegrator) {
203 this.starter = starterIntegrator;
204 }
205
206 /** Start the integration.
207 * <p>This method computes one step using the underlying starter integrator,
208 * and initializes the Nordsieck vector at step start. The starter integrator
209 * purpose is only to establish initial conditions, it does not really change
210 * time by itself. The top level multistep integrator remains in charge of
211 * handling time propagation and events handling as it will starts its own
212 * computation right from the beginning. In a sense, the starter integrator
213 * can be seen as a dummy one and so it will never trigger any user event nor
214 * call any user step handler.</p>
215 * @param equations complete set of differential equations to integrate
216 * @param initialState initial state (time, primary and secondary state vectors)
217 * @param t target time for the integration
218 * (can be set to a value smaller than <code>t0</code> for backward integration)
219 * @exception DimensionMismatchException if arrays dimension do not match equations settings
220 * @exception NumberIsTooSmallException if integration step is too small
221 * @exception MaxCountExceededException if the number of functions evaluations is exceeded
222 * @exception NoBracketingException if the location of an event cannot be bracketed
223 */
224 protected void start(final FieldExpandableODE<T> equations, final FieldODEState<T> initialState, final T t)
225 throws DimensionMismatchException, NumberIsTooSmallException,
226 MaxCountExceededException, NoBracketingException {
227
228 // make sure NO user event nor user step handler is triggered,
229 // this is the task of the top level integrator, not the task
230 // of the starter integrator
231 starter.clearEventHandlers();
232 starter.clearStepHandlers();
233
234 // set up one specific step handler to extract initial Nordsieck vector
235 starter.addStepHandler(new FieldNordsieckInitializer(equations.getMapper(), (nSteps + 3) / 2));
236
237 // start integration, expecting a InitializationCompletedMarkerException
238 try {
239
240 starter.integrate(equations, initialState, t);
241
242 // we should not reach this step
243 throw new MathIllegalStateException(LocalizedFormats.MULTISTEP_STARTER_STOPPED_EARLY);
244
245 } catch (InitializationCompletedMarkerException icme) { // NOPMD
246 // this is the expected nominal interruption of the start integrator
247
248 // count the evaluations used by the starter
249 getEvaluationsCounter().increment(starter.getEvaluations());
250
251 }
252
253 // remove the specific step handler
254 starter.clearStepHandlers();
255
256 }
257
258 /** Initialize the high order scaled derivatives at step start.
259 * @param h step size to use for scaling
260 * @param t first steps times
261 * @param y first steps states
262 * @param yDot first steps derivatives
263 * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
264 * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
265 */
266 protected abstract Array2DRowFieldMatrix<T> initializeHighOrderDerivatives(final T h, final T[] t,
267 final T[][] y,
268 final T[][] yDot);
269
270 /** Get the minimal reduction factor for stepsize control.
271 * @return minimal reduction factor
272 */
273 public double getMinReduction() {
274 return minReduction;
275 }
276
277 /** Set the minimal reduction factor for stepsize control.
278 * @param minReduction minimal reduction factor
279 */
280 public void setMinReduction(final double minReduction) {
281 this.minReduction = minReduction;
282 }
283
284 /** Get the maximal growth factor for stepsize control.
285 * @return maximal growth factor
286 */
287 public double getMaxGrowth() {
288 return maxGrowth;
289 }
290
291 /** Set the maximal growth factor for stepsize control.
292 * @param maxGrowth maximal growth factor
293 */
294 public void setMaxGrowth(final double maxGrowth) {
295 this.maxGrowth = maxGrowth;
296 }
297
298 /** Get the safety factor for stepsize control.
299 * @return safety factor
300 */
301 public double getSafety() {
302 return safety;
303 }
304
305 /** Set the safety factor for stepsize control.
306 * @param safety safety factor
307 */
308 public void setSafety(final double safety) {
309 this.safety = safety;
310 }
311
312 /** Get the number of steps of the multistep method (excluding the one being computed).
313 * @return number of steps of the multistep method (excluding the one being computed)
314 */
315 public int getNSteps() {
316 return nSteps;
317 }
318
319 /** Rescale the instance.
320 * <p>Since the scaled and Nordsieck arrays are shared with the caller,
321 * this method has the side effect of rescaling this arrays in the caller too.</p>
322 * @param newStepSize new step size to use in the scaled and Nordsieck arrays
323 */
324 protected void rescale(final T newStepSize) {
325
326 final T ratio = newStepSize.divide(getStepSize());
327 for (int i = 0; i < scaled.length; ++i) {
328 scaled[i] = scaled[i].multiply(ratio);
329 }
330
331 final T[][] nData = nordsieck.getDataRef();
332 T power = ratio;
333 for (int i = 0; i < nData.length; ++i) {
334 power = power.multiply(ratio);
335 final T[] nDataI = nData[i];
336 for (int j = 0; j < nDataI.length; ++j) {
337 nDataI[j] = nDataI[j].multiply(power);
338 }
339 }
340
341 setStepSize(newStepSize);
342
343 }
344
345
346 /** Compute step grow/shrink factor according to normalized error.
347 * @param error normalized error of the current step
348 * @return grow/shrink factor for next step
349 */
350 protected T computeStepGrowShrinkFactor(final T error) {
351 return MathUtils.min(error.getField().getZero().add(maxGrowth),
352 MathUtils.max(error.getField().getZero().add(minReduction),
353 error.pow(exp).multiply(safety)));
354 }
355
356 /** Specialized step handler storing the first step.
357 */
358 private class FieldNordsieckInitializer implements FieldStepHandler<T> {
359
360 /** Equation mapper. */
361 private final FieldEquationsMapper<T> mapper;
362
363 /** Steps counter. */
364 private int count;
365
366 /** Saved start. */
367 private FieldODEStateAndDerivative<T> savedStart;
368
369 /** First steps times. */
370 private final T[] t;
371
372 /** First steps states. */
373 private final T[][] y;
374
375 /** First steps derivatives. */
376 private final T[][] yDot;
377
378 /** Simple constructor.
379 * @param mapper equation mapper
380 * @param nbStartPoints number of start points (including the initial point)
381 */
382 FieldNordsieckInitializer(final FieldEquationsMapper<T> mapper, final int nbStartPoints) {
383 this.mapper = mapper;
384 this.count = 0;
385 this.t = MathArrays.buildArray(getField(), nbStartPoints);
386 this.y = MathArrays.buildArray(getField(), nbStartPoints, -1);
387 this.yDot = MathArrays.buildArray(getField(), nbStartPoints, -1);
388 }
389
390 /** {@inheritDoc} */
391 public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast)
392 throws MaxCountExceededException {
393
394
395 if (count == 0) {
396 // first step, we need to store also the point at the beginning of the step
397 final FieldODEStateAndDerivative<T> prev = interpolator.getPreviousState();
398 savedStart = prev;
399 t[count] = prev.getTime();
400 y[count] = mapper.mapState(prev);
401 yDot[count] = mapper.mapDerivative(prev);
402 }
403
404 // store the point at the end of the step
405 ++count;
406 final FieldODEStateAndDerivative<T> curr = interpolator.getCurrentState();
407 t[count] = curr.getTime();
408 y[count] = mapper.mapState(curr);
409 yDot[count] = mapper.mapDerivative(curr);
410
411 if (count == t.length - 1) {
412
413 // this was the last point we needed, we can compute the derivatives
414 setStepSize(t[t.length - 1].subtract(t[0]).divide(t.length - 1));
415
416 // first scaled derivative
417 scaled = MathArrays.buildArray(getField(), yDot[0].length);
418 for (int j = 0; j < scaled.length; ++j) {
419 scaled[j] = yDot[0][j].multiply(getStepSize());
420 }
421
422 // higher order derivatives
423 nordsieck = initializeHighOrderDerivatives(getStepSize(), t, y, yDot);
424
425 // stop the integrator now that all needed steps have been handled
426 setStepStart(savedStart);
427 throw new InitializationCompletedMarkerException();
428
429 }
430
431 }
432
433 /** {@inheritDoc} */
434 public void init(final FieldODEStateAndDerivative<T> initialState, T finalTime) {
435 // nothing to do
436 }
437
438 }
439
440 /** Marker exception used ONLY to stop the starter integrator after first step. */
441 private static class InitializationCompletedMarkerException
442 extends RuntimeException {
443
444 /** Serializable version identifier. */
445 private static final long serialVersionUID = -1914085471038046418L;
446
447 /** Simple constructor. */
448 InitializationCompletedMarkerException() {
449 super((Throwable) null);
450 }
451
452 }
453
454}
Note: See TracBrowser for help on using the repository browser.