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.nonstiff.interpolators.RungeKuttaStateInterpolator;
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 FixedStepRungeKuttaIntegrator extends AbstractIntegrator
55          implements ExplicitRungeKuttaIntegrator {
56  
57      /** Time steps from Butcher array (without the first zero). */
58      private final double[] c;
59  
60      /** Internal weights from Butcher array (without the first empty row). */
61      private final double[][] a;
62  
63      /** External weights for the high order method from Butcher array. */
64      private final double[] b;
65  
66      /** Integration step. */
67      private final double step;
68  
69      /** Simple constructor.
70       * Build a Runge-Kutta integrator with the given
71       * step. The default step handler does nothing.
72       * @param name name of the method
73       * @param step integration step
74       */
75      protected FixedStepRungeKuttaIntegrator(final String name, final double step) {
76          super(name);
77          this.c    = getC();
78          this.a    = getA();
79          this.b    = getB();
80          this.step = FastMath.abs(step);
81      }
82  
83      /** Getter for the default, positive step-size assigned at constructor level.
84       * @return step
85       */
86      public double getDefaultStep() {
87          return this.step;
88      }
89  
90      /** Create an interpolator.
91       * @param forward integration direction indicator
92       * @param yDotK slopes at the intermediate points
93       * @param globalPreviousState start of the global step
94       * @param globalCurrentState end of the global step
95       * @param mapper equations mapper for the all equations
96       * @return external weights for the high order method from Butcher array
97       */
98      protected abstract RungeKuttaStateInterpolator createInterpolator(boolean forward, double[][] yDotK,
99                                                                        ODEStateAndDerivative globalPreviousState,
100                                                                       ODEStateAndDerivative globalCurrentState,
101                                                                       EquationsMapper mapper);
102 
103     /** {@inheritDoc} */
104     @Override
105     public ODEStateAndDerivative integrate(final ExpandableODE equations,
106                                            final ODEState initialState, final double finalTime)
107         throws MathIllegalArgumentException, MathIllegalStateException {
108 
109         sanityChecks(initialState, finalTime);
110         setStepStart(initIntegration(equations, initialState, finalTime));
111         final boolean forward = finalTime > initialState.getTime();
112 
113         // create some internal working arrays
114         final int        stages = c.length + 1;
115         double[]         y      = getStepStart().getCompleteState();
116         final double[][] yDotK  = new double[stages][];
117         final double[]   yTmp   = new double[y.length];
118 
119         // set up integration control objects
120         if (forward) {
121             if (getStepStart().getTime() + step >= finalTime) {
122                 setStepSize(finalTime - getStepStart().getTime());
123             } else {
124                 setStepSize(step);
125             }
126         } else {
127             if (getStepStart().getTime() - step <= finalTime) {
128                 setStepSize(finalTime - getStepStart().getTime());
129             } else {
130                 setStepSize(-step);
131             }
132         }
133 
134         // main integration loop
135         setIsLastStep(false);
136         do {
137 
138             // first stage
139             y        = getStepStart().getCompleteState();
140             yDotK[0] = getStepStart().getCompleteDerivative();
141 
142             // next stages
143             ExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(), getStepStart().getTime(), y,
144                     getStepSize(), a, c, yDotK);
145 
146             incrementEvaluations(stages - 1);
147 
148             // estimate the state at the end of the step
149             for (int j = 0; j < y.length; ++j) {
150                 double sum    = b[0] * yDotK[0][j];
151                 for (int l = 1; l < stages; ++l) {
152                     sum    += b[l] * yDotK[l][j];
153                 }
154                 yTmp[j] = y[j] + getStepSize() * sum;
155                 if (Double.isNaN(yTmp[j])) {
156                     throw new MathIllegalStateException(LocalizedODEFormats.NAN_APPEARING_DURING_INTEGRATION,
157                                                         getStepStart().getTime() + getStepSize());
158                 }
159 
160             }
161             final double stepEnd   = getStepStart().getTime() + getStepSize();
162             final double[] yDotTmp = computeDerivatives(stepEnd, yTmp);
163             final ODEStateAndDerivative stateTmp =
164                 equations.getMapper().mapStateAndDerivative(stepEnd, yTmp, yDotTmp);
165 
166             // discrete events handling
167             System.arraycopy(yTmp, 0, y, 0, y.length);
168             setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp,
169                                                        equations.getMapper()),
170                                     finalTime));
171 
172             if (!isLastStep()) {
173 
174                 // stepsize control for next step
175                 final double  nextT      = getStepStart().getTime() + getStepSize();
176                 final boolean nextIsLast = forward ? (nextT >= finalTime) : (nextT <= finalTime);
177                 if (nextIsLast) {
178                     setStepSize(finalTime - getStepStart().getTime());
179                 }
180             }
181 
182         } while (!isLastStep());
183 
184         final ODEStateAndDerivative finalState = getStepStart();
185         setStepStart(null);
186         setStepSize(Double.NaN);
187         return finalState;
188 
189     }
190 
191 }