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