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