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