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.math3.ode.nonstiff;
19 
20 import org.apache.commons.math3.Field;
21 import org.apache.commons.math3.RealFieldElement;
22 import org.apache.commons.math3.ode.FieldEquationsMapper;
23 import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
24 
25 /**
26  * This class represents an interpolator over the last step during an
27  * ODE integration for the 6th order Luther integrator.
28  *
29  * <p>This interpolator computes dense output inside the last
30  * step computed. The interpolation equation is consistent with the
31  * integration scheme.</p>
32  *
33  * @see LutherFieldIntegrator
34  * @param <T> the type of the field elements
35  * @since 3.6
36  */
37 
38 class LutherFieldStepInterpolator<T extends RealFieldElement<T>>
39     extends RungeKuttaFieldStepInterpolator<T> {
40 
41     /** -49 - 49 q. */
42     private final T c5a;
43 
44     /** 392 + 287 q. */
45     private final T c5b;
46 
47     /** -637 - 357 q. */
48     private final T c5c;
49 
50     /** 833 + 343 q. */
51     private final T c5d;
52 
53     /** -49 + 49 q. */
54     private final T c6a;
55 
56     /** -392 - 287 q. */
57     private final T c6b;
58 
59     /** -637 + 357 q. */
60     private final T c6c;
61 
62     /** 833 - 343 q. */
63     private final T c6d;
64 
65     /** 49 + 49 q. */
66     private final T d5a;
67 
68     /** -1372 - 847 q. */
69     private final T d5b;
70 
71     /** 2254 + 1029 q */
72     private final T d5c;
73 
74     /** 49 - 49 q. */
75     private final T d6a;
76 
77     /** -1372 + 847 q. */
78     private final T d6b;
79 
80     /** 2254 - 1029 q */
81     private final T d6c;
82 
83     /** Simple constructor.
84      * @param field field to which the time and state vector elements belong
85      * @param forward integration direction indicator
86      * @param yDotK slopes at the intermediate points
87      * @param globalPreviousState start of the global step
88      * @param globalCurrentState end of the global step
89      * @param softPreviousState start of the restricted step
90      * @param softCurrentState end of the restricted step
91      * @param mapper equations mapper for the all equations
92      */
LutherFieldStepInterpolator(final Field<T> field, final boolean forward, final T[][] yDotK, final FieldODEStateAndDerivative<T> globalPreviousState, final FieldODEStateAndDerivative<T> globalCurrentState, final FieldODEStateAndDerivative<T> softPreviousState, final FieldODEStateAndDerivative<T> softCurrentState, final FieldEquationsMapper<T> mapper)93     LutherFieldStepInterpolator(final Field<T> field, final boolean forward,
94                                 final T[][] yDotK,
95                                 final FieldODEStateAndDerivative<T> globalPreviousState,
96                                 final FieldODEStateAndDerivative<T> globalCurrentState,
97                                 final FieldODEStateAndDerivative<T> softPreviousState,
98                                 final FieldODEStateAndDerivative<T> softCurrentState,
99                                 final FieldEquationsMapper<T> mapper) {
100         super(field, forward, yDotK,
101               globalPreviousState, globalCurrentState, softPreviousState, softCurrentState,
102               mapper);
103         final T q = field.getZero().add(21).sqrt();
104         c5a = q.multiply(  -49).add(  -49);
105         c5b = q.multiply(  287).add(  392);
106         c5c = q.multiply( -357).add( -637);
107         c5d = q.multiply(  343).add(  833);
108         c6a = q.multiply(   49).add(  -49);
109         c6b = q.multiply( -287).add(  392);
110         c6c = q.multiply(  357).add( -637);
111         c6d = q.multiply( -343).add(  833);
112         d5a = q.multiply(   49).add(   49);
113         d5b = q.multiply( -847).add(-1372);
114         d5c = q.multiply( 1029).add( 2254);
115         d6a = q.multiply(  -49).add(   49);
116         d6b = q.multiply(  847).add(-1372);
117         d6c = q.multiply(-1029).add( 2254);
118     }
119 
120     /** {@inheritDoc} */
121     @Override
create(final Field<T> newField, final boolean newForward, final T[][] newYDotK, final FieldODEStateAndDerivative<T> newGlobalPreviousState, final FieldODEStateAndDerivative<T> newGlobalCurrentState, final FieldODEStateAndDerivative<T> newSoftPreviousState, final FieldODEStateAndDerivative<T> newSoftCurrentState, final FieldEquationsMapper<T> newMapper)122     protected LutherFieldStepInterpolator<T> create(final Field<T> newField, final boolean newForward, final T[][] newYDotK,
123                                                     final FieldODEStateAndDerivative<T> newGlobalPreviousState,
124                                                     final FieldODEStateAndDerivative<T> newGlobalCurrentState,
125                                                     final FieldODEStateAndDerivative<T> newSoftPreviousState,
126                                                     final FieldODEStateAndDerivative<T> newSoftCurrentState,
127                                                     final FieldEquationsMapper<T> newMapper) {
128         return new LutherFieldStepInterpolator<T>(newField, newForward, newYDotK,
129                                                   newGlobalPreviousState, newGlobalCurrentState,
130                                                   newSoftPreviousState, newSoftCurrentState,
131                                                   newMapper);
132     }
133 
134     /** {@inheritDoc} */
135     @SuppressWarnings("unchecked")
136     @Override
computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper, final T time, final T theta, final T thetaH, final T oneMinusThetaH)137     protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
138                                                                                    final T time, final T theta,
139                                                                                    final T thetaH, final T oneMinusThetaH) {
140 
141         // the coefficients below have been computed by solving the
142         // order conditions from a theorem from Butcher (1963), using
143         // the method explained in Folkmar Bornemann paper "Runge-Kutta
144         // Methods, Trees, and Maple", Center of Mathematical Sciences, Munich
145         // University of Technology, February 9, 2001
146         //<http://wwwzenger.informatik.tu-muenchen.de/selcuk/sjam012101.html>
147 
148         // the method is implemented in the rkcheck tool
149         // <https://www.spaceroots.org/software/rkcheck/index.html>.
150         // Running it for order 5 gives the following order conditions
151         // for an interpolator:
152         // order 1 conditions
153         // \sum_{i=1}^{i=s}\left(b_{i} \right) =1
154         // order 2 conditions
155         // \sum_{i=1}^{i=s}\left(b_{i} c_{i}\right) = \frac{\theta}{2}
156         // order 3 conditions
157         // \sum_{i=2}^{i=s}\left(b_{i} \sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j} \right)}\right) = \frac{\theta^{2}}{6}
158         // \sum_{i=1}^{i=s}\left(b_{i} c_{i}^{2}\right) = \frac{\theta^{2}}{3}
159         // order 4 conditions
160         // \sum_{i=3}^{i=s}\left(b_{i} \sum_{j=2}^{j=i-1}{\left(a_{i,j} \sum_{k=1}^{k=j-1}{\left(a_{j,k} c_{k} \right)} \right)}\right) = \frac{\theta^{3}}{24}
161         // \sum_{i=2}^{i=s}\left(b_{i} \sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j}^{2} \right)}\right) = \frac{\theta^{3}}{12}
162         // \sum_{i=2}^{i=s}\left(b_{i} c_{i}\sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j} \right)}\right) = \frac{\theta^{3}}{8}
163         // \sum_{i=1}^{i=s}\left(b_{i} c_{i}^{3}\right) = \frac{\theta^{3}}{4}
164         // order 5 conditions
165         // \sum_{i=4}^{i=s}\left(b_{i} \sum_{j=3}^{j=i-1}{\left(a_{i,j} \sum_{k=2}^{k=j-1}{\left(a_{j,k} \sum_{l=1}^{l=k-1}{\left(a_{k,l} c_{l} \right)} \right)} \right)}\right) = \frac{\theta^{4}}{120}
166         // \sum_{i=3}^{i=s}\left(b_{i} \sum_{j=2}^{j=i-1}{\left(a_{i,j} \sum_{k=1}^{k=j-1}{\left(a_{j,k} c_{k}^{2} \right)} \right)}\right) = \frac{\theta^{4}}{60}
167         // \sum_{i=3}^{i=s}\left(b_{i} \sum_{j=2}^{j=i-1}{\left(a_{i,j} c_{j}\sum_{k=1}^{k=j-1}{\left(a_{j,k} c_{k} \right)} \right)}\right) = \frac{\theta^{4}}{40}
168         // \sum_{i=2}^{i=s}\left(b_{i} \sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j}^{3} \right)}\right) = \frac{\theta^{4}}{20}
169         // \sum_{i=3}^{i=s}\left(b_{i} c_{i}\sum_{j=2}^{j=i-1}{\left(a_{i,j} \sum_{k=1}^{k=j-1}{\left(a_{j,k} c_{k} \right)} \right)}\right) = \frac{\theta^{4}}{30}
170         // \sum_{i=2}^{i=s}\left(b_{i} c_{i}\sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j}^{2} \right)}\right) = \frac{\theta^{4}}{15}
171         // \sum_{i=2}^{i=s}\left(b_{i} \left(\sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j} \right)} \right)^{2}\right) = \frac{\theta^{4}}{20}
172         // \sum_{i=2}^{i=s}\left(b_{i} c_{i}^{2}\sum_{j=1}^{j=i-1}{\left(a_{i,j} c_{j} \right)}\right) = \frac{\theta^{4}}{10}
173         // \sum_{i=1}^{i=s}\left(b_{i} c_{i}^{4}\right) = \frac{\theta^{4}}{5}
174 
175         // The a_{j,k} and c_{k} are given by the integrator Butcher arrays. What remains to solve
176         // are the b_i for the interpolator. They are found by solving the above equations.
177         // For a given interpolator, some equations are redundant, so in our case when we select
178         // all equations from order 1 to 4, we still don't have enough independent equations
179         // to solve from b_1 to b_7. We need to also select one equation from order 5. Here,
180         // we selected the last equation. It appears this choice implied at least the last 3 equations
181         // are fulfilled, but some of the former ones are not, so the resulting interpolator is order 5.
182         // At the end, we get the b_i as polynomials in theta.
183 
184         final T coeffDot1 =  theta.multiply(theta.multiply(theta.multiply(theta.multiply(   21        ).add( -47          )).add(   36         )).add( -54     /   5.0)).add(1);
185         final T coeffDot2 =  time.getField().getZero();
186         final T coeffDot3 =  theta.multiply(theta.multiply(theta.multiply(theta.multiply(  112        ).add(-608    /  3.0)).add(  320   / 3.0 )).add(-208    /  15.0));
187         final T coeffDot4 =  theta.multiply(theta.multiply(theta.multiply(theta.multiply( -567  /  5.0).add( 972    /  5.0)).add( -486   / 5.0 )).add( 324    /  25.0));
188         final T coeffDot5 =  theta.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(5)).add(c5b.divide(15))).add(c5c.divide(30))).add(c5d.divide(150)));
189         final T coeffDot6 =  theta.multiply(theta.multiply(theta.multiply(theta.multiply(c6a.divide(5)).add(c6b.divide(15))).add(c6c.divide(30))).add(c6d.divide(150)));
190         final T coeffDot7 =  theta.multiply(theta.multiply(theta.multiply(                                             3.0 ).add(   -3         )).add(   3   /   5.0));
191         final T[] interpolatedState;
192         final T[] interpolatedDerivatives;
193 
194         if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
195 
196             final T s         = thetaH;
197             final T coeff1    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(  21    /  5.0).add( -47    /  4.0)).add(   12         )).add( -27    /   5.0)).add(1));
198             final T coeff2    = time.getField().getZero();
199             final T coeff3    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply( 112    /  5.0).add(-152    /  3.0)).add(  320   / 9.0 )).add(-104    /  15.0)));
200             final T coeff4    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(-567    / 25.0).add( 243    /  5.0)).add( -162   / 5.0 )).add( 162    /  25.0)));
201             final T coeff5    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(25)).add(c5b.divide(60))).add(c5c.divide(90))).add(c5d.divide(300))));
202             final T coeff6    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(c6a.divide(25)).add(c6b.divide(60))).add(c6c.divide(90))).add(c6d.divide(300))));
203             final T coeff7    = s.multiply(theta.multiply(theta.multiply(theta.multiply(                                      3    /  4.0 ).add(   -1         )).add(   3    /  10.0)));
204             interpolatedState       = previousStateLinearCombination(coeff1, coeff2, coeff3, coeff4, coeff5, coeff6, coeff7);
205             interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4, coeffDot5, coeffDot6, coeffDot7);
206         } else {
207 
208             final T s         = oneMinusThetaH;
209             final T coeff1    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply( -21   /   5.0).add(   151  /  20.0)).add( -89   /   20.0)).add(  19 /  20.0)).add(- 1 / 20.0));
210             final T coeff2    = time.getField().getZero();
211             final T coeff3    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(-112   /   5.0).add(   424  /  15.0)).add( -328  /   45.0)).add( -16 /  45.0)).add(-16 /  45.0));
212             final T coeff4    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply( 567   /  25.0).add(  -648  /  25.0)).add(  162  /   25.0))));
213             final T coeff5    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(d5a.divide(25)).add(d5b.divide(300))).add(d5c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0));
214             final T coeff6    = s.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(d6a.divide(25)).add(d6b.divide(300))).add(d6c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0));
215             final T coeff7    = s.multiply(               theta.multiply(theta.multiply(theta.multiply(                        -3  /   4.0 ).add(   1   /    4.0)).add(  -1 /  20.0)).add( -1 /  20.0));
216             interpolatedState       = currentStateLinearCombination(coeff1, coeff2, coeff3, coeff4, coeff5, coeff6, coeff7);
217             interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4, coeffDot5, coeffDot6, coeffDot7);
218         }
219 
220         return new FieldODEStateAndDerivative<T>(time, interpolatedState, interpolatedDerivatives);
221 
222     }
223 
224 }
225