View Javadoc
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 }