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 agents.anac.y2019.harddealer.math3.ode;
|
---|
19 |
|
---|
20 | import agents.anac.y2019.harddealer.math3.Field;
|
---|
21 | import agents.anac.y2019.harddealer.math3.RealFieldElement;
|
---|
22 | import agents.anac.y2019.harddealer.math3.exception.DimensionMismatchException;
|
---|
23 | import agents.anac.y2019.harddealer.math3.exception.MathIllegalStateException;
|
---|
24 | import agents.anac.y2019.harddealer.math3.exception.MaxCountExceededException;
|
---|
25 | import agents.anac.y2019.harddealer.math3.exception.NoBracketingException;
|
---|
26 | import agents.anac.y2019.harddealer.math3.exception.NumberIsTooSmallException;
|
---|
27 | import agents.anac.y2019.harddealer.math3.exception.util.LocalizedFormats;
|
---|
28 | import agents.anac.y2019.harddealer.math3.linear.Array2DRowFieldMatrix;
|
---|
29 | import agents.anac.y2019.harddealer.math3.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
|
---|
30 | import agents.anac.y2019.harddealer.math3.ode.nonstiff.DormandPrince853FieldIntegrator;
|
---|
31 | import agents.anac.y2019.harddealer.math3.ode.sampling.FieldStepHandler;
|
---|
32 | import agents.anac.y2019.harddealer.math3.ode.sampling.FieldStepInterpolator;
|
---|
33 | import agents.anac.y2019.harddealer.math3.util.FastMath;
|
---|
34 | import agents.anac.y2019.harddealer.math3.util.MathArrays;
|
---|
35 | import 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 | */
|
---|
68 | public 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 | }
|
---|