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  import org.hipparchus.CalculusFieldElement;
26  import org.hipparchus.Field;
27  import org.hipparchus.ode.FieldEquationsMapper;
28  import org.hipparchus.ode.FieldODEStateAndDerivative;
29  
30  /**
31   * This class represents an interpolator over the last step during an
32   * ODE integration for the 5(4) Dormand-Prince integrator.
33   *
34   * @see DormandPrince54Integrator
35   *
36   * @param <T> the type of the field elements
37   */
38  
39  class DormandPrince54FieldStateInterpolator<T extends CalculusFieldElement<T>>
40      extends RungeKuttaFieldStateInterpolator<T> {
41  
42      /** Last row of the Butcher-array internal weights, element 0. */
43      private final T a70;
44  
45      // element 1 is zero, so it is neither stored nor used
46  
47      /** Last row of the Butcher-array internal weights, element 2. */
48      private final T a72;
49  
50      /** Last row of the Butcher-array internal weights, element 3. */
51      private final T a73;
52  
53      /** Last row of the Butcher-array internal weights, element 4. */
54      private final T a74;
55  
56      /** Last row of the Butcher-array internal weights, element 5. */
57      private final T a75;
58  
59      /** Shampine (1986) Dense output, element 0. */
60      private final T d0;
61  
62      // element 1 is zero, so it is neither stored nor used
63  
64      /** Shampine (1986) Dense output, element 2. */
65      private final T d2;
66  
67      /** Shampine (1986) Dense output, element 3. */
68      private final T d3;
69  
70      /** Shampine (1986) Dense output, element 4. */
71      private final T d4;
72  
73      /** Shampine (1986) Dense output, element 5. */
74      private final T d5;
75  
76      /** Shampine (1986) Dense output, element 6. */
77      private final T d6;
78  
79      /** Simple constructor.
80       * @param field field to which the time and state vector elements belong
81       * @param forward integration direction indicator
82       * @param yDotK slopes at the intermediate points
83       * @param globalPreviousState start of the global step
84       * @param globalCurrentState end of the global step
85       * @param softPreviousState start of the restricted step
86       * @param softCurrentState end of the restricted step
87       * @param mapper equations mapper for the all equations
88       */
89      DormandPrince54FieldStateInterpolator(final Field<T> field, final boolean forward,
90                                            final T[][] yDotK,
91                                            final FieldODEStateAndDerivative<T> globalPreviousState,
92                                            final FieldODEStateAndDerivative<T> globalCurrentState,
93                                            final FieldODEStateAndDerivative<T> softPreviousState,
94                                            final FieldODEStateAndDerivative<T> softCurrentState,
95                                            final FieldEquationsMapper<T> mapper) {
96          super(field, forward, yDotK,
97                globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
98                mapper);
99          final T one = field.getOne();
100         a70 = one.newInstance(   35.0 /  384.0);
101         a72 = one.newInstance(  500.0 / 1113.0);
102         a73 = one.newInstance(  125.0 /  192.0);
103         a74 = one.newInstance(-2187.0 / 6784.0);
104         a75 = one.newInstance(   11.0 /   84.0);
105         d0  = one.newInstance(-12715105075.0 /  11282082432.0);
106         d2  = one.newInstance( 87487479700.0 /  32700410799.0);
107         d3  = one.newInstance(-10690763975.0 /   1880347072.0);
108         d4  = one.newInstance(701980252875.0 / 199316789632.0);
109         d5  = one.newInstance( -1453857185.0 /    822651844.0);
110         d6  = one.newInstance(    69997945.0 /     29380423.0);
111     }
112 
113     /** {@inheritDoc} */
114     @Override
115     protected DormandPrince54FieldStateInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK,
116                                                               final FieldODEStateAndDerivative<T> newGlobalPreviousState,
117                                                               final FieldODEStateAndDerivative<T> newGlobalCurrentState,
118                                                               final FieldODEStateAndDerivative<T> newSoftPreviousState,
119                                                               final FieldODEStateAndDerivative<T> newSoftCurrentState,
120                                                               final FieldEquationsMapper<T> newMapper) {
121         return new DormandPrince54FieldStateInterpolator<T>(newField, newForward, newYDotK,
122                                                             newGlobalPreviousState, newGlobalCurrentState,
123                                                             newSoftPreviousState, newSoftCurrentState,
124                                                             newMapper);
125     }
126     /** {@inheritDoc} */
127     @SuppressWarnings("unchecked")
128     @Override
129     protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
130                                                                                    final T time, final T theta,
131                                                                                    final T thetaH, final T oneMinusThetaH) {
132 
133         // interpolate
134         final T one      = time.getField().getOne();
135         final T eta      = one.subtract(theta);
136         final T twoTheta = theta.multiply(2);
137         final T dot2     = one.subtract(twoTheta);
138         final T dot3     = theta.multiply(theta.multiply(-3).add(2));
139         final T dot4     = twoTheta.multiply(theta.multiply(twoTheta.subtract(3)).add(1));
140         final T[] interpolatedState;
141         final T[] interpolatedDerivatives;
142         if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
143             final T f1        = thetaH;
144             final T f2        = f1.multiply(eta);
145             final T f3        = f2.multiply(theta);
146             final T f4        = f3.multiply(eta);
147             final T coeff0    = f1.multiply(a70).
148                                 subtract(f2.multiply(a70.subtract(1))).
149                                 add(f3.multiply(a70.multiply(2).subtract(1))).
150                                 add(f4.multiply(d0));
151             final T coeff1    = time.getField().getZero();
152             final T coeff2    = f1.multiply(a72).
153                                 subtract(f2.multiply(a72)).
154                                 add(f3.multiply(a72.multiply(2))).
155                                 add(f4.multiply(d2));
156             final T coeff3    = f1.multiply(a73).
157                                 subtract(f2.multiply(a73)).
158                                 add(f3.multiply(a73.multiply(2))).
159                                 add(f4.multiply(d3));
160             final T coeff4    = f1.multiply(a74).
161                                 subtract(f2.multiply(a74)).
162                                 add(f3.multiply(a74.multiply(2))).
163                                 add(f4.multiply(d4));
164             final T coeff5    = f1.multiply(a75).
165                                 subtract(f2.multiply(a75)).
166                                 add(f3.multiply(a75.multiply(2))).
167                                 add(f4.multiply(d5));
168             final T coeff6    = f4.multiply(d6).subtract(f3);
169             final T coeffDot0 = a70.
170                                 subtract(dot2.multiply(a70.subtract(1))).
171                                 add(dot3.multiply(a70.multiply(2).subtract(1))).
172                                 add(dot4.multiply(d0));
173             final T coeffDot1 = time.getField().getZero();
174             final T coeffDot2 = a72.
175                                 subtract(dot2.multiply(a72)).
176                                 add(dot3.multiply(a72.multiply(2))).
177                                 add(dot4.multiply(d2));
178             final T coeffDot3 = a73.
179                                 subtract(dot2.multiply(a73)).
180                                 add(dot3.multiply(a73.multiply(2))).
181                                 add(dot4.multiply(d3));
182             final T coeffDot4 = a74.
183                                 subtract(dot2.multiply(a74)).
184                                 add(dot3.multiply(a74.multiply(2))).
185                                 add(dot4.multiply(d4));
186             final T coeffDot5 = a75.
187                                 subtract(dot2.multiply(a75)).
188                                 add(dot3.multiply(a75.multiply(2))).
189                                 add(dot4.multiply(d5));
190             final T coeffDot6 = dot4.multiply(d6).subtract(dot3);
191             interpolatedState       = previousStateLinearCombination(coeff0, coeff1, coeff2, coeff3,
192                                                                      coeff4, coeff5, coeff6);
193             interpolatedDerivatives = derivativeLinearCombination(coeffDot0, coeffDot1, coeffDot2, coeffDot3,
194                                                                   coeffDot4, coeffDot5, coeffDot6);
195         } else {
196             final T f1        = oneMinusThetaH.negate();
197             final T f2        = oneMinusThetaH.multiply(theta);
198             final T f3        = f2.multiply(theta);
199             final T f4        = f3.multiply(eta);
200             final T coeff0    = f1.multiply(a70).
201                                 subtract(f2.multiply(a70.subtract(1))).
202                                 add(f3.multiply(a70.multiply(2).subtract(1))).
203                                 add(f4.multiply(d0));
204             final T coeff1    = time.getField().getZero();
205             final T coeff2    = f1.multiply(a72).
206                                 subtract(f2.multiply(a72)).
207                                 add(f3.multiply(a72.multiply(2))).
208                                 add(f4.multiply(d2));
209             final T coeff3    = f1.multiply(a73).
210                                 subtract(f2.multiply(a73)).
211                                 add(f3.multiply(a73.multiply(2))).
212                                 add(f4.multiply(d3));
213             final T coeff4    = f1.multiply(a74).
214                                 subtract(f2.multiply(a74)).
215                                 add(f3.multiply(a74.multiply(2))).
216                                 add(f4.multiply(d4));
217             final T coeff5    = f1.multiply(a75).
218                                 subtract(f2.multiply(a75)).
219                                 add(f3.multiply(a75.multiply(2))).
220                                 add(f4.multiply(d5));
221             final T coeff6    = f4.multiply(d6).subtract(f3);
222             final T coeffDot0 = a70.
223                                 subtract(dot2.multiply(a70.subtract(1))).
224                                 add(dot3.multiply(a70.multiply(2).subtract(1))).
225                                 add(dot4.multiply(d0));
226             final T coeffDot1 = time.getField().getZero();
227             final T coeffDot2 = a72.
228                                 subtract(dot2.multiply(a72)).
229                                 add(dot3.multiply(a72.multiply(2))).
230                                 add(dot4.multiply(d2));
231             final T coeffDot3 = a73.
232                                 subtract(dot2.multiply(a73)).
233                                 add(dot3.multiply(a73.multiply(2))).
234                                 add(dot4.multiply(d3));
235             final T coeffDot4 = a74.
236                                 subtract(dot2.multiply(a74)).
237                                 add(dot3.multiply(a74.multiply(2))).
238                                 add(dot4.multiply(d4));
239             final T coeffDot5 = a75.
240                                 subtract(dot2.multiply(a75)).
241                                 add(dot3.multiply(a75.multiply(2))).
242                                 add(dot4.multiply(d5));
243             final T coeffDot6 = dot4.multiply(d6).subtract(dot3);
244             interpolatedState       = currentStateLinearCombination(coeff0, coeff1, coeff2, coeff3,
245                                                                     coeff4, coeff5, coeff6);
246             interpolatedDerivatives = derivativeLinearCombination(coeffDot0, coeffDot1, coeffDot2, coeffDot3,
247                                                                   coeffDot4, coeffDot5, coeffDot6);
248         }
249         return mapper.mapStateAndDerivative(time, interpolatedState, interpolatedDerivatives);
250 
251     }
252 
253 }