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