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.nonstiff.interpolators.RungeKuttaFieldStateInterpolator;
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 FixedStepRungeKuttaFieldIntegrator<T extends CalculusFieldElement<T>>
62      extends AbstractFieldIntegrator<T> implements FieldExplicitRungeKuttaIntegrator<T> {
63  
64      /** Time steps from Butcher array (without the first zero). */
65      private final T[] c;
66  
67      /** Internal weights from Butcher array (without the first empty row). */
68      private final T[][] a;
69  
70      /** External weights for the high order method from Butcher array. */
71      private final T[] b;
72  
73      /** Time steps from Butcher array (without the first zero). */
74      private double[] realC = new double[0];
75  
76      /** Internal weights from Butcher array (without the first empty row). Real version, optional. */
77      private double[][] realA = new double[0][];
78  
79      /** External weights for the high order method from Butcher array. Real version, optional. */
80      private double[] realB = new double[0];
81  
82      /** Integration step. */
83      private final T step;
84  
85      /** Flag setting whether coefficients in Butcher array are interpreted as Field or real numbers. */
86      private boolean usingFieldCoefficients;
87  
88      /** Simple constructor.
89       * Build a Runge-Kutta integrator with the given
90       * step. The default step handler does nothing.
91       * @param field field to which the time and state vector elements belong
92       * @param name name of the method
93       * @param step integration step
94       */
95      protected FixedStepRungeKuttaFieldIntegrator(final Field<T> field, final String name, final T step) {
96          super(field, name);
97          this.c    = getC();
98          this.a    = getA();
99          this.b    = getB();
100         this.step = step.abs();
101         this.usingFieldCoefficients = false;
102     }
103 
104     /** Getter for the default, positive step-size assigned at constructor level.
105      * @return step
106      */
107     public T getDefaultStep() {
108         return this.step;
109     }
110 
111     /**
112      * Setter for the flag between real or Field coefficients in the Butcher array.
113      *
114      * @param usingFieldCoefficients new value for flag
115      */
116     public void setUsingFieldCoefficients(boolean usingFieldCoefficients) {
117         this.usingFieldCoefficients = usingFieldCoefficients;
118     }
119 
120     /** {@inheritDoc} */
121     @Override
122     public boolean isUsingFieldCoefficients() {
123         return usingFieldCoefficients;
124     }
125 
126     /** {@inheritDoc} */
127     @Override
128     public int getNumberOfStages() {
129         return b.length;
130     }
131 
132     /** Create an interpolator.
133      * @param forward integration direction indicator
134      * @param yDotK slopes at the intermediate points
135      * @param globalPreviousState start of the global step
136      * @param globalCurrentState end of the global step
137      * @param mapper equations mapper for the all equations
138      * @return external weights for the high order method from Butcher array
139      */
140     protected abstract RungeKuttaFieldStateInterpolator<T> createInterpolator(boolean forward, T[][] yDotK,
141                                                                               FieldODEStateAndDerivative<T> globalPreviousState,
142                                                                               FieldODEStateAndDerivative<T> globalCurrentState,
143                                                                               FieldEquationsMapper<T> mapper);
144 
145     /** {@inheritDoc} */
146     @Override
147     protected FieldODEStateAndDerivative<T> initIntegration(FieldExpandableODE<T> eqn, FieldODEState<T> s0, T t) {
148         if (!isUsingFieldCoefficients()) {
149             realA = getRealA();
150             realB = getRealB();
151             realC = getRealC();
152         }
153         return super.initIntegration(eqn, s0, t);
154     }
155 
156     /** {@inheritDoc} */
157     @Override
158     public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations,
159                                                    final FieldODEState<T> initialState, final T finalTime)
160         throws MathIllegalArgumentException, MathIllegalStateException {
161 
162         sanityChecks(initialState, finalTime);
163         setStepStart(initIntegration(equations, initialState, finalTime));
164         final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0;
165 
166         // create some internal working arrays
167         final int   stages = getNumberOfStages();
168         final T[][] yDotK  = MathArrays.buildArray(getField(), stages, -1);
169         MathArrays.buildArray(getField(), equations.getMapper().getTotalDimension());
170 
171         // set up integration control objects
172         if (forward) {
173             if (getStepStart().getTime().add(step).subtract(finalTime).getReal() >= 0) {
174                 setStepSize(finalTime.subtract(getStepStart().getTime()));
175             } else {
176                 setStepSize(step);
177             }
178         } else {
179             if (getStepStart().getTime().subtract(step).subtract(finalTime).getReal() <= 0) {
180                 setStepSize(finalTime.subtract(getStepStart().getTime()));
181             } else {
182                 setStepSize(step.negate());
183             }
184         }
185 
186         // main integration loop
187         setIsLastStep(false);
188         do {
189 
190             // first stage
191             final T[] y = getStepStart().getCompleteState();
192             yDotK[0]    = getStepStart().getCompleteDerivative();
193 
194             // next stages
195             final T[] yTmp;
196             if (isUsingFieldCoefficients()) {
197                 FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(), getStepStart().getTime(),
198                         y, getStepSize(), a, c, yDotK);
199                 yTmp = FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y, yDotK, getStepSize(), b);
200             } else {
201                 FieldExplicitRungeKuttaIntegrator.applyInternalButcherWeights(getEquations(), getStepStart().getTime(),
202                         y, getStepSize(), realA, realC, yDotK);
203                 yTmp = FieldExplicitRungeKuttaIntegrator.applyExternalButcherWeights(y, yDotK, getStepSize(), realB);
204             }
205 
206             incrementEvaluations(stages - 1);
207 
208             final T stepEnd   = getStepStart().getTime().add(getStepSize());
209             final T[] yDotTmp = computeDerivatives(stepEnd, yTmp);
210             final FieldODEStateAndDerivative<T> stateTmp = equations.getMapper().mapStateAndDerivative(stepEnd, yTmp, yDotTmp);
211 
212             // discrete events handling
213             setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp, equations.getMapper()),
214                                     finalTime));
215 
216             if (!isLastStep()) {
217 
218                 // stepsize control for next step
219                 final T  nextT      = getStepStart().getTime().add(getStepSize());
220                 final boolean nextIsLast = forward ?
221                                            (nextT.subtract(finalTime).getReal() >= 0) :
222                                            (nextT.subtract(finalTime).getReal() <= 0);
223                 if (nextIsLast) {
224                     setStepSize(finalTime.subtract(getStepStart().getTime()));
225                 }
226             }
227 
228         } while (!isLastStep());
229 
230         final FieldODEStateAndDerivative<T> finalState = getStepStart();
231         setStepStart(null);
232         setStepSize(null);
233         return finalState;
234 
235     }
236 
237 }