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 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 */
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
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<>(newField, newForward, newYDotK,
129 newGlobalPreviousState, newGlobalCurrentState,
130 newSoftPreviousState, newSoftCurrentState,
131 newMapper);
132 }
133
134 /** {@inheritDoc} */
135 @SuppressWarnings("unchecked")
136 @Override
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<>(time, interpolatedState, interpolatedDerivatives);
221 }
222 }