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