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