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 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 */ DormandPrince54FieldStepInterpolator(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)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 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)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<T>(newField, newForward, newYDotK, 118 newGlobalPreviousState, newGlobalCurrentState, 119 newSoftPreviousState, newSoftCurrentState, 120 newMapper); 121 } 122 /** {@inheritDoc} */ 123 @SuppressWarnings("unchecked") 124 @Override computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper, final T time, final T theta, final T thetaH, final T oneMinusThetaH)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<T>(time, interpolatedState, interpolatedDerivatives); 246 247 } 248 249 } 250