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 * https://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 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22
23 package org.hipparchus.ode.nonstiff;
24
25 import org.hipparchus.Field;
26 import org.hipparchus.CalculusFieldElement;
27 import org.hipparchus.exception.MathIllegalArgumentException;
28 import org.hipparchus.exception.MathIllegalStateException;
29 import org.hipparchus.linear.Array2DRowFieldMatrix;
30 import org.hipparchus.linear.FieldMatrix;
31 import org.hipparchus.ode.FieldEquationsMapper;
32 import org.hipparchus.ode.FieldExpandableODE;
33 import org.hipparchus.ode.FieldODEState;
34 import org.hipparchus.ode.FieldODEStateAndDerivative;
35 import org.hipparchus.ode.LocalizedODEFormats;
36 import org.hipparchus.ode.MultistepFieldIntegrator;
37 import org.hipparchus.util.MathArrays;
38
39
40 /** Base class for {@link AdamsBashforthFieldIntegrator Adams-Bashforth} and
41 * {@link AdamsMoultonFieldIntegrator Adams-Moulton} integrators.
42 * @param <T> the type of the field elements
43 */
44 public abstract class AdamsFieldIntegrator<T extends CalculusFieldElement<T>> extends MultistepFieldIntegrator<T> {
45
46 /** Transformer. */
47 private final AdamsNordsieckFieldTransformer<T> transformer;
48
49 /**
50 * Build an Adams integrator with the given order and step control parameters.
51 * @param field field to which the time and state vector elements belong
52 * @param name name of the method
53 * @param nSteps number of steps of the method excluding the one being computed
54 * @param order order of the method
55 * @param minStep minimal step (sign is irrelevant, regardless of
56 * integration direction, forward or backward), the last step can
57 * be smaller than this
58 * @param maxStep maximal step (sign is irrelevant, regardless of
59 * integration direction, forward or backward), the last step can
60 * be smaller than this
61 * @param scalAbsoluteTolerance allowed absolute error
62 * @param scalRelativeTolerance allowed relative error
63 * @exception MathIllegalArgumentException if order is 1 or less
64 */
65 public AdamsFieldIntegrator(final Field<T> field, final String name,
66 final int nSteps, final int order,
67 final double minStep, final double maxStep,
68 final double scalAbsoluteTolerance,
69 final double scalRelativeTolerance)
70 throws MathIllegalArgumentException {
71 super(field, name, nSteps, order, minStep, maxStep,
72 scalAbsoluteTolerance, scalRelativeTolerance);
73 transformer = AdamsNordsieckFieldTransformer.getInstance(field, nSteps);
74 }
75
76 /**
77 * Build an Adams integrator with the given order and step control parameters.
78 * @param field field to which the time and state vector elements belong
79 * @param name name of the method
80 * @param nSteps number of steps of the method excluding the one being computed
81 * @param order order of the method
82 * @param minStep minimal step (sign is irrelevant, regardless of
83 * integration direction, forward or backward), the last step can
84 * be smaller than this
85 * @param maxStep maximal step (sign is irrelevant, regardless of
86 * integration direction, forward or backward), the last step can
87 * be smaller than this
88 * @param vecAbsoluteTolerance allowed absolute error
89 * @param vecRelativeTolerance allowed relative error
90 * @exception IllegalArgumentException if order is 1 or less
91 */
92 public AdamsFieldIntegrator(final Field<T> field, final String name,
93 final int nSteps, final int order,
94 final double minStep, final double maxStep,
95 final double[] vecAbsoluteTolerance,
96 final double[] vecRelativeTolerance)
97 throws IllegalArgumentException {
98 super(field, name, nSteps, order, minStep, maxStep,
99 vecAbsoluteTolerance, vecRelativeTolerance);
100 transformer = AdamsNordsieckFieldTransformer.getInstance(field, nSteps);
101 }
102
103 /** {@inheritDoc} */
104 @Override
105 public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
106 final FieldODEState<T> initialState,
107 final T finalTime)
108 throws MathIllegalArgumentException, MathIllegalStateException {
109
110 sanityChecks(initialState, finalTime);
111 setStepStart(initIntegration(equations, initialState, finalTime));
112 final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0;
113
114 // compute the initial Nordsieck vector using the configured starter integrator
115 start(equations, getStepStart(), finalTime);
116
117 // reuse the step that was chosen by the starter integrator
118 FieldODEStateAndDerivative<T> stepStart = getStepStart();
119 FieldODEStateAndDerivative<T> stepEnd =
120 AdamsFieldStateInterpolator.taylor(equations.getMapper(), stepStart,
121 stepStart.getTime().add(getStepSize()),
122 getStepSize(), scaled, nordsieck);
123
124 // main integration loop
125 setIsLastStep(false);
126 final T[] y = stepStart.getCompleteState();
127 do {
128
129 T[] predictedY = null;
130 final T[] predictedScaled = MathArrays.buildArray(getField(), y.length);
131 Array2DRowFieldMatrix<T> predictedNordsieck = null;
132 double error = 10;
133 while (error >= 1.0) {
134
135 // predict a first estimate of the state at step end
136 predictedY = stepEnd.getCompleteState();
137
138 // evaluate the derivative
139 final T[] yDot = computeDerivatives(stepEnd.getTime(), predictedY);
140
141 // predict Nordsieck vector at step end
142 for (int j = 0; j < predictedScaled.length; ++j) {
143 predictedScaled[j] = getStepSize().multiply(yDot[j]);
144 }
145 predictedNordsieck = updateHighOrderDerivativesPhase1(nordsieck);
146 updateHighOrderDerivativesPhase2(scaled, predictedScaled, predictedNordsieck);
147
148 // evaluate error
149 error = errorEstimation(y, stepEnd.getTime(), predictedY, predictedScaled, predictedNordsieck);
150 if (Double.isNaN(error)) {
151 throw new MathIllegalStateException(LocalizedODEFormats.NAN_APPEARING_DURING_INTEGRATION,
152 stepEnd.getTime().getReal());
153 }
154
155 if (error >= 1.0) {
156 // reject the step and attempt to reduce error by stepsize control
157 final double factor = computeStepGrowShrinkFactor(error);
158 rescale(getStepSizeHelper().filterStep(getStepSize().multiply(factor), forward, false));
159 stepEnd = AdamsFieldStateInterpolator.taylor(equations.getMapper(), getStepStart(),
160 getStepStart().getTime().add(getStepSize()),
161 getStepSize(),
162 scaled,
163 nordsieck);
164
165 }
166 }
167
168 final AdamsFieldStateInterpolator<T> interpolator =
169 finalizeStep(getStepSize(), predictedY, predictedScaled, predictedNordsieck,
170 forward, getStepStart(), stepEnd, equations.getMapper());
171
172 // discrete events handling
173 setStepStart(acceptStep(interpolator, finalTime));
174 scaled = interpolator.getScaled();
175 nordsieck = interpolator.getNordsieck();
176
177 if (!isLastStep()) {
178
179 if (resetOccurred()) {
180
181 // some events handler has triggered changes that
182 // invalidate the derivatives, we need to restart from scratch
183 start(equations, getStepStart(), finalTime);
184
185 final T nextT = getStepStart().getTime().add(getStepSize());
186 final boolean nextIsLast = forward ?
187 nextT.subtract(finalTime).getReal() >= 0 :
188 nextT.subtract(finalTime).getReal() <= 0;
189 final T hNew = nextIsLast ? finalTime.subtract(getStepStart().getTime()) : getStepSize();
190
191 rescale(hNew);
192 System.arraycopy(getStepStart().getCompleteState(), 0, y, 0, y.length);
193
194 } else {
195
196 // stepsize control for next step
197 final double factor = computeStepGrowShrinkFactor(error);
198 final T scaledH = getStepSize().multiply(factor);
199 final T nextT = getStepStart().getTime().add(scaledH);
200 final boolean nextIsLast = forward ?
201 nextT.subtract(finalTime).getReal() >= 0 :
202 nextT.subtract(finalTime).getReal() <= 0;
203 T hNew = getStepSizeHelper().filterStep(scaledH, forward, nextIsLast);
204
205 final T filteredNextT = getStepStart().getTime().add(hNew);
206 final boolean filteredNextIsLast = forward ?
207 filteredNextT.subtract(finalTime).getReal() >= 0 :
208 filteredNextT.subtract(finalTime).getReal() <= 0;
209 if (filteredNextIsLast) {
210 hNew = finalTime.subtract(getStepStart().getTime());
211 }
212
213 rescale(hNew);
214 System.arraycopy(predictedY, 0, y, 0, y.length);
215
216 }
217
218 stepEnd = AdamsFieldStateInterpolator.taylor(equations.getMapper(), getStepStart(),
219 getStepStart().getTime().add(getStepSize()),
220 getStepSize(), scaled, nordsieck);
221
222 }
223
224 } while (!isLastStep());
225
226 final FieldODEStateAndDerivative<T> finalState = getStepStart();
227 setStepStart(null);
228 setStepSize(null);
229 return finalState;
230
231 }
232
233 /** {@inheritDoc} */
234 @Override
235 protected Array2DRowFieldMatrix<T> initializeHighOrderDerivatives(final T h, final T[] t,
236 final T[][] y,
237 final T[][] yDot) {
238 return transformer.initializeHighOrderDerivatives(h, t, y, yDot);
239 }
240
241 /** Update the high order scaled derivatives for Adams integrators (phase 1).
242 * <p>The complete update of high order derivatives has a form similar to:
243 * \[
244 * r_{n+1} = (s_1(n) - s_1(n+1)) P^{-1} u + P^{-1} A P r_n
245 * \]
246 * this method computes the P<sup>-1</sup> A P r<sub>n</sub> part.</p>
247 * @param highOrder high order scaled derivatives
248 * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
249 * @return updated high order derivatives
250 * @see #updateHighOrderDerivativesPhase2(CalculusFieldElement[], CalculusFieldElement[], Array2DRowFieldMatrix)
251 */
252 public Array2DRowFieldMatrix<T> updateHighOrderDerivativesPhase1(final Array2DRowFieldMatrix<T> highOrder) {
253 return transformer.updateHighOrderDerivativesPhase1(highOrder);
254 }
255
256 /** Update the high order scaled derivatives Adams integrators (phase 2).
257 * <p>The complete update of high order derivatives has a form similar to:
258 * \[
259 * r_{n+1} = (s_1(n) - s_1(n+1)) P^{-1} u + P^{-1} A P r_n
260 * \]
261 * this method computes the (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u part.</p>
262 * <p>Phase 1 of the update must already have been performed.</p>
263 * @param start first order scaled derivatives at step start
264 * @param end first order scaled derivatives at step end
265 * @param highOrder high order scaled derivatives, will be modified
266 * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
267 * @see #updateHighOrderDerivativesPhase1(Array2DRowFieldMatrix)
268 */
269 public void updateHighOrderDerivativesPhase2(final T[] start, final T[] end,
270 final Array2DRowFieldMatrix<T> highOrder) {
271 transformer.updateHighOrderDerivativesPhase2(start, end, highOrder);
272 }
273
274 /** Estimate error.
275 * @param previousState state vector at step start
276 * @param predictedTime time at step end
277 * @param predictedState predicted state vector at step end
278 * @param predictedScaled predicted value of the scaled derivatives at step end
279 * @param predictedNordsieck predicted value of the Nordsieck vector at step end
280 * @return estimated normalized local discretization error
281 * @since 2.0
282 */
283 protected abstract double errorEstimation(T[] previousState, T predictedTime,
284 T[] predictedState, T[] predictedScaled,
285 FieldMatrix<T> predictedNordsieck);
286
287 /** Finalize the step.
288 * @param stepSize step size used in the scaled and Nordsieck arrays
289 * @param predictedState predicted state at end of step
290 * @param predictedScaled predicted first scaled derivative
291 * @param predictedNordsieck predicted Nordsieck vector
292 * @param isForward integration direction indicator
293 * @param globalPreviousState start of the global step
294 * @param globalCurrentState end of the global step
295 * @param equationsMapper mapper for ODE equations primary and secondary components
296 * @return step interpolator
297 * @since 2.0
298 */
299 protected abstract AdamsFieldStateInterpolator<T> finalizeStep(T stepSize, T[] predictedState,
300 T[] predictedScaled, Array2DRowFieldMatrix<T> predictedNordsieck,
301 boolean isForward,
302 FieldODEStateAndDerivative<T> globalPreviousState,
303 FieldODEStateAndDerivative<T> globalCurrentState,
304 FieldEquationsMapper<T> equationsMapper);
305
306 }