1 /*
2 * Licensed to the Hipparchus project 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 Hipparchus project 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 package org.hipparchus.ode.nonstiff.interpolators;
19
20 import org.hipparchus.ode.EquationsMapper;
21 import org.hipparchus.ode.ODEStateAndDerivative;
22 import org.hipparchus.ode.nonstiff.LutherIntegrator;
23 import org.hipparchus.util.FastMath;
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 LutherIntegrator
34 */
35
36 public class LutherStateInterpolator extends RungeKuttaStateInterpolator {
37
38 /** Serializable version identifier */
39 private static final long serialVersionUID = 20160328;
40
41 /** Square root. */
42 private static final double Q = FastMath.sqrt(21);
43
44 /** Simple constructor.
45 * @param forward integration direction indicator
46 * @param yDotK slopes at the intermediate points
47 * @param globalPreviousState start of the global step
48 * @param globalCurrentState end of the global step
49 * @param softPreviousState start of the restricted step
50 * @param softCurrentState end of the restricted step
51 * @param mapper equations mapper for the all equations
52 */
53 public LutherStateInterpolator(final boolean forward,
54 final double[][] yDotK,
55 final ODEStateAndDerivative globalPreviousState,
56 final ODEStateAndDerivative globalCurrentState,
57 final ODEStateAndDerivative softPreviousState,
58 final ODEStateAndDerivative softCurrentState,
59 final EquationsMapper mapper) {
60 super(forward, yDotK, globalPreviousState, globalCurrentState, softPreviousState, softCurrentState, mapper);
61 }
62
63 /** {@inheritDoc} */
64 @Override
65 protected LutherStateInterpolator create(final boolean newForward, final double[][] newYDotK,
66 final ODEStateAndDerivative newGlobalPreviousState,
67 final ODEStateAndDerivative newGlobalCurrentState,
68 final ODEStateAndDerivative newSoftPreviousState,
69 final ODEStateAndDerivative newSoftCurrentState,
70 final EquationsMapper newMapper) {
71 return new LutherStateInterpolator(newForward, newYDotK,
72 newGlobalPreviousState, newGlobalCurrentState,
73 newSoftPreviousState, newSoftCurrentState,
74 newMapper);
75 }
76
77 /** {@inheritDoc} */
78 @Override
79 protected ODEStateAndDerivative computeInterpolatedStateAndDerivatives(final EquationsMapper mapper,
80 final double time, final double theta,
81 final double thetaH, final double oneMinusThetaH) {
82
83 // the coefficients below have been computed by solving the
84 // order conditions from a theorem from Butcher (1963), using
85 // the method explained in Folkmar Bornemann paper "Runge-Kutta
86 // Methods, Trees, and Maple", Center of Mathematical Sciences, Munich
87 // University of Technology, February 9, 2001
88 //<http://wwwzenger.informatik.tu-muenchen.de/selcuk/sjam012101.html>
89
90 // the method is implemented in the rkcheck tool
91 // <https://www.spaceroots.org/software/rkcheck/index.html>.
92 // Running it for order 5 gives the following order conditions
93 // for an interpolator:
94 // order 1 conditions
95 // \sum_{i=1}^{i=s}\left(b_{i} \right) =1
96 // order 2 conditions
97 // \sum_{i=1}^{i=s}\left(b_{i} c_{i}\right) = \frac{\theta}{2}
98 // order 3 conditions
99 // \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}
100 // \sum_{i=1}^{i=s}\left(b_{i} c_{i}^{2}\right) = \frac{\theta^{2}}{3}
101 // order 4 conditions
102 // \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}
103 // \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}
104 // \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}
105 // \sum_{i=1}^{i=s}\left(b_{i} c_{i}^{3}\right) = \frac{\theta^{3}}{4}
106 // order 5 conditions
107 // \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}
108 // \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}
109 // \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}
110 // \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}
111 // \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}
112 // \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}
113 // \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}
114 // \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}
115 // \sum_{i=1}^{i=s}\left(b_{i} c_{i}^{4}\right) = \frac{\theta^{4}}{5}
116
117 // The a_{j,k} and c_{k} are given by the integrator Butcher arrays. What remains to solve
118 // are the b_i for the interpolator. They are found by solving the above equations.
119 // For a given interpolator, some equations are redundant, so in our case when we select
120 // all equations from order 1 to 4, we still don't have enough independent equations
121 // to solve from b_1 to b_7. We need to also select one equation from order 5. Here,
122 // we selected the last equation. It appears this choice implied at least the last 3 equations
123 // are fulfilled, but some of the former ones are not, so the resulting interpolator is order 5.
124 // At the end, we get the b_i as polynomials in theta.
125
126 final double[] interpolatedState;
127 final double[] interpolatedDerivatives;
128
129 final double coeffDot1 = 1 + theta * ( -54 / 5.0 + theta * ( 36 + theta * ( -47 + theta * 21)));
130 final double coeffDot2 = 0;
131 final double coeffDot3 = theta * (-208 / 15.0 + theta * ( 320 / 3.0 + theta * (-608 / 3.0 + theta * 112)));
132 final double coeffDot4 = theta * ( 324 / 25.0 + theta * ( -486 / 5.0 + theta * ( 972 / 5.0 + theta * -567 / 5.0)));
133 final double coeffDot5 = theta * ((833 + 343 * Q) / 150.0 + theta * ((-637 - 357 * Q) / 30.0 + theta * ((392 + 287 * Q) / 15.0 + theta * (-49 - 49 * Q) / 5.0)));
134 final double coeffDot6 = theta * ((833 - 343 * Q) / 150.0 + theta * ((-637 + 357 * Q) / 30.0 + theta * ((392 - 287 * Q) / 15.0 + theta * (-49 + 49 * Q) / 5.0)));
135 final double coeffDot7 = theta * ( 3 / 5.0 + theta * ( -3 + theta * 3));
136
137 if (getGlobalPreviousState() != null && theta <= 0.5) {
138
139 final double coeff1 = 1 + theta * ( -27 / 5.0 + theta * ( 12 + theta * ( -47 / 4.0 + theta * 21 / 5.0)));
140 final double coeff2 = 0;
141 final double coeff3 = theta * (-104 / 15.0 + theta * ( 320 / 9.0 + theta * (-152 / 3.0 + theta * 112 / 5.0)));
142 final double coeff4 = theta * ( 162 / 25.0 + theta * ( -162 / 5.0 + theta * ( 243 / 5.0 + theta * -567 / 25.0)));
143 final double coeff5 = theta * ((833 + 343 * Q) / 300.0 + theta * ((-637 - 357 * Q) / 90.0 + theta * ((392 + 287 * Q) / 60.0 + theta * (-49 - 49 * Q) / 25.0)));
144 final double coeff6 = theta * ((833 - 343 * Q) / 300.0 + theta * ((-637 + 357 * Q) / 90.0 + theta * ((392 - 287 * Q) / 60.0 + theta * (-49 + 49 * Q) / 25.0)));
145 final double coeff7 = theta * ( 3 / 10.0 + theta * ( -1 + theta * ( 3 / 4.0)));
146 interpolatedState = previousStateLinearCombination(thetaH * coeff1, thetaH * coeff2,
147 thetaH * coeff3, thetaH * coeff4,
148 thetaH * coeff5, thetaH * coeff6,
149 thetaH * coeff7);
150 interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4, coeffDot5, coeffDot6, coeffDot7);
151 } else {
152
153 final double coeff1 = -1 / 20.0 + theta * ( 19 / 20.0 + theta * ( -89 / 20.0 + theta * ( 151 / 20.0 + theta * -21 / 5.0)));
154 final double coeff2 = 0;
155 final double coeff3 = -16 / 45.0 + theta * ( -16 / 45.0 + theta * ( -328 / 45.0 + theta * ( 424 / 15.0 + theta * -112 / 5.0)));
156 final double coeff4 = theta * ( theta * ( 162 / 25.0 + theta * ( -648 / 25.0 + theta * 567 / 25.0)));
157 final double coeff5 = -49 / 180.0 + theta * ( -49 / 180.0 + theta * ((2254 + 1029 * Q) / 900.0 + theta * ((-1372 - 847 * Q) / 300.0 + theta * ( 49 + 49 * Q) / 25.0)));
158 final double coeff6 = -49 / 180.0 + theta * ( -49 / 180.0 + theta * ((2254 - 1029 * Q) / 900.0 + theta * ((-1372 + 847 * Q) / 300.0 + theta * ( 49 - 49 * Q) / 25.0)));
159 final double coeff7 = -1 / 20.0 + theta * ( -1 / 20.0 + theta * ( 1 / 4.0 + theta * ( -3 / 4.0)));
160 interpolatedState = currentStateLinearCombination(oneMinusThetaH * coeff1, oneMinusThetaH * coeff2,
161 oneMinusThetaH * coeff3, oneMinusThetaH * coeff4,
162 oneMinusThetaH * coeff5, oneMinusThetaH * coeff6,
163 oneMinusThetaH * coeff7);
164 interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4, coeffDot5, coeffDot6, coeffDot7);
165 }
166
167 return mapper.mapStateAndDerivative(time, interpolatedState, interpolatedDerivatives);
168
169 }
170
171 }