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  
26  import org.hipparchus.CalculusFieldElement;
27  import org.hipparchus.Field;
28  import org.hipparchus.exception.MathIllegalArgumentException;
29  import org.hipparchus.exception.MathIllegalStateException;
30  import org.hipparchus.ode.AbstractFieldIntegrator;
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.FieldOrdinaryDifferentialEquation;
36  import org.hipparchus.util.MathArrays;
37  
38  /**
39   * This class implements the common part of all fixed step Runge-Kutta
40   * integrators for Ordinary Differential Equations.
41   *
42   * <p>These methods are explicit Runge-Kutta methods, their Butcher
43   * arrays are as follows :</p>
44   * <pre>
45   *    0  |
46   *   c2  | a21
47   *   c3  | a31  a32
48   *   ... |        ...
49   *   cs  | as1  as2  ...  ass-1
50   *       |--------------------------
51   *       |  b1   b2  ...   bs-1  bs
52   * </pre>
53   *
54   * @see EulerFieldIntegrator
55   * @see ClassicalRungeKuttaFieldIntegrator
56   * @see GillFieldIntegrator
57   * @see MidpointFieldIntegrator
58   * @param <T> the type of the field elements
59   */
60  
61  public abstract class RungeKuttaFieldIntegrator<T extends CalculusFieldElement<T>>
62      extends AbstractFieldIntegrator<T>
63      implements FieldButcherArrayProvider<T> {
64  
65      /** Time steps from Butcher array (without the first zero). */
66      private final T[] c;
67  
68      /** Internal weights from Butcher array (without the first empty row). */
69      private final T[][] a;
70  
71      /** External weights for the high order method from Butcher array. */
72      private final T[] b;
73  
74      /** Integration step. */
75      private final T step;
76  
77      /** Simple constructor.
78       * Build a Runge-Kutta integrator with the given
79       * step. The default step handler does nothing.
80       * @param field field to which the time and state vector elements belong
81       * @param name name of the method
82       * @param step integration step
83       */
84      protected RungeKuttaFieldIntegrator(final Field<T> field, final String name, final T step) {
85          super(field, name);
86          this.c    = getC();
87          this.a    = getA();
88          this.b    = getB();
89          this.step = step.abs();
90      }
91  
92      /** Getter for the default, positive step-size assigned at constructor level.
93       * @return step
94       */
95      public T getDefaultStep() {
96          return this.step;
97      }
98  
99      /** Create a fraction.
100      * @param p numerator
101      * @param q denominator
102      * @return p/q computed in the instance field
103      */
104     protected T fraction(final int p, final int q) {
105         return getField().getZero().add(p).divide(q);
106     }
107 
108     /** Create an interpolator.
109      * @param forward integration direction indicator
110      * @param yDotK slopes at the intermediate points
111      * @param globalPreviousState start of the global step
112      * @param globalCurrentState end of the global step
113      * @param mapper equations mapper for the all equations
114      * @return external weights for the high order method from Butcher array
115      */
116     protected abstract RungeKuttaFieldStateInterpolator<T> createInterpolator(boolean forward, T[][] yDotK,
117                                                                              FieldODEStateAndDerivative<T> globalPreviousState,
118                                                                              FieldODEStateAndDerivative<T> globalCurrentState,
119                                                                              FieldEquationsMapper<T> mapper);
120 
121     /** {@inheritDoc} */
122     @Override
123     public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
124                                                    final FieldODEState<T> initialState, final T finalTime)
125         throws MathIllegalArgumentException, MathIllegalStateException {
126 
127         sanityChecks(initialState, finalTime);
128         setStepStart(initIntegration(equations, initialState, finalTime));
129         final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0;
130 
131         // create some internal working arrays
132         final int   stages = c.length + 1;
133         final T[][] yDotK  = MathArrays.buildArray(getField(), stages, -1);
134         final T[]   yTmp   = MathArrays.buildArray(getField(), equations.getMapper().getTotalDimension());
135 
136         // set up integration control objects
137         if (forward) {
138             if (getStepStart().getTime().add(step).subtract(finalTime).getReal() >= 0) {
139                 setStepSize(finalTime.subtract(getStepStart().getTime()));
140             } else {
141                 setStepSize(step);
142             }
143         } else {
144             if (getStepStart().getTime().subtract(step).subtract(finalTime).getReal() <= 0) {
145                 setStepSize(finalTime.subtract(getStepStart().getTime()));
146             } else {
147                 setStepSize(step.negate());
148             }
149         }
150 
151         // main integration loop
152         setIsLastStep(false);
153         do {
154 
155             // first stage
156             final T[] y = getStepStart().getCompleteState();
157             yDotK[0]    = getStepStart().getCompleteDerivative();
158 
159             // next stages
160             for (int k = 1; k < stages; ++k) {
161 
162                 for (int j = 0; j < y.length; ++j) {
163                     T sum = yDotK[0][j].multiply(a[k-1][0]);
164                     for (int l = 1; l < k; ++l) {
165                         sum = sum.add(yDotK[l][j].multiply(a[k-1][l]));
166                     }
167                     yTmp[j] = y[j].add(getStepSize().multiply(sum));
168                 }
169 
170                 yDotK[k] = computeDerivatives(getStepStart().getTime().add(getStepSize().multiply(c[k-1])), yTmp);
171 
172             }
173 
174             // estimate the state at the end of the step
175             for (int j = 0; j < y.length; ++j) {
176                 T sum = yDotK[0][j].multiply(b[0]);
177                 for (int l = 1; l < stages; ++l) {
178                     sum = sum.add(yDotK[l][j].multiply(b[l]));
179                 }
180                 yTmp[j] = y[j].add(getStepSize().multiply(sum));
181             }
182             final T stepEnd   = getStepStart().getTime().add(getStepSize());
183             final T[] yDotTmp = computeDerivatives(stepEnd, yTmp);
184             final FieldODEStateAndDerivative<T> stateTmp = equations.getMapper().mapStateAndDerivative(stepEnd, yTmp, yDotTmp);
185 
186             // discrete events handling
187             setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp, equations.getMapper()),
188                                     finalTime));
189 
190             if (!isLastStep()) {
191 
192                 // stepsize control for next step
193                 final T  nextT      = getStepStart().getTime().add(getStepSize());
194                 final boolean nextIsLast = forward ?
195                                            (nextT.subtract(finalTime).getReal() >= 0) :
196                                            (nextT.subtract(finalTime).getReal() <= 0);
197                 if (nextIsLast) {
198                     setStepSize(finalTime.subtract(getStepStart().getTime()));
199                 }
200             }
201 
202         } while (!isLastStep());
203 
204         final FieldODEStateAndDerivative<T> finalState = getStepStart();
205         setStepStart(null);
206         setStepSize(null);
207         return finalState;
208 
209     }
210 
211     /** Fast computation of a single step of ODE integration.
212      * <p>This method is intended for the limited use case of
213      * very fast computation of only one step without using any of the
214      * rich features of general integrators that may take some time
215      * to set up (i.e. no step handlers, no events handlers, no additional
216      * states, no interpolators, no error control, no evaluations count,
217      * no sanity checks ...). It handles the strict minimum of computation,
218      * so it can be embedded in outer loops.</p>
219      * <p>
220      * This method is <em>not</em> used at all by the {@link #integrate(FieldExpandableODE,
221      * FieldODEState, CalculusFieldElement)} method. It also completely ignores the step set at
222      * construction time, and uses only a single step to go from {@code t0} to {@code t}.
223      * </p>
224      * <p>
225      * As this method does not use any of the state-dependent features of the integrator,
226      * it should be reasonably thread-safe <em>if and only if</em> the provided differential
227      * equations are themselves thread-safe.
228      * </p>
229      * @param equations differential equations to integrate
230      * @param t0 initial time
231      * @param y0 initial value of the state vector at t0
232      * @param t target time for the integration
233      * (can be set to a value smaller than {@code t0} for backward integration)
234      * @return state vector at {@code t}
235      */
236     public T[] singleStep(final FieldOrdinaryDifferentialEquation<T> equations,
237                           final T t0, final T[] y0, final T t) {
238 
239         // create some internal working arrays
240         final T[] y       = y0.clone();
241         final int stages  = c.length + 1;
242         final T[][] yDotK = MathArrays.buildArray(getField(), stages, -1);
243         final T[] yTmp    = y0.clone();
244 
245         // first stage
246         final T h = t.subtract(t0);
247         yDotK[0] = equations.computeDerivatives(t0, y);
248 
249         // next stages
250         for (int k = 1; k < stages; ++k) {
251 
252             for (int j = 0; j < y0.length; ++j) {
253                 T sum = yDotK[0][j].multiply(a[k-1][0]);
254                 for (int l = 1; l < k; ++l) {
255                     sum = sum.add(yDotK[l][j].multiply(a[k-1][l]));
256                 }
257                 yTmp[j] = y[j].add(h.multiply(sum));
258             }
259 
260             yDotK[k] = equations.computeDerivatives(t0.add(h.multiply(c[k-1])), yTmp);
261 
262         }
263 
264         // estimate the state at the end of the step
265         for (int j = 0; j < y0.length; ++j) {
266             T sum = yDotK[0][j].multiply(b[0]);
267             for (int l = 1; l < stages; ++l) {
268                 sum = sum.add(yDotK[l][j].multiply(b[l]));
269             }
270             y[j] = y[j].add(h.multiply(sum));
271         }
272 
273         return y;
274 
275     }
276 
277 }