001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018package org.apache.commons.math4.legacy.ode.nonstiff; 019 020import org.apache.commons.math4.legacy.core.Field; 021import org.apache.commons.math4.legacy.core.RealFieldElement; 022import org.apache.commons.math4.legacy.exception.DimensionMismatchException; 023import org.apache.commons.math4.legacy.exception.MaxCountExceededException; 024import org.apache.commons.math4.legacy.exception.NoBracketingException; 025import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException; 026import org.apache.commons.math4.legacy.ode.FieldEquationsMapper; 027import org.apache.commons.math4.legacy.ode.FieldExpandableODE; 028import org.apache.commons.math4.legacy.ode.FieldODEState; 029import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative; 030import org.apache.commons.math4.legacy.core.MathArrays; 031 032/** 033 * This class implements the common part of all embedded Runge-Kutta 034 * integrators for Ordinary Differential Equations. 035 * 036 * <p>These methods are embedded explicit Runge-Kutta methods with two 037 * sets of coefficients allowing to estimate the error, their Butcher 038 * arrays are as follows : 039 * <pre> 040 * 0 | 041 * c2 | a21 042 * c3 | a31 a32 043 * ... | ... 044 * cs | as1 as2 ... ass-1 045 * |-------------------------- 046 * | b1 b2 ... bs-1 bs 047 * | b'1 b'2 ... b's-1 b's 048 * </pre> 049 * 050 * <p>In fact, we rather use the array defined by ej = bj - b'j to 051 * compute directly the error rather than computing two estimates and 052 * then comparing them.</p> 053 * 054 * <p>Some methods are qualified as <i>fsal</i> (first same as last) 055 * methods. This means the last evaluation of the derivatives in one 056 * step is the same as the first in the next step. Then, this 057 * evaluation can be reused from one step to the next one and the cost 058 * of such a method is really s-1 evaluations despite the method still 059 * has s stages. This behaviour is true only for successful steps, if 060 * the step is rejected after the error estimation phase, no 061 * evaluation is saved. For an <i>fsal</i> method, we have cs = 1 and 062 * asi = bi for all i.</p> 063 * 064 * @param <T> the type of the field elements 065 * @since 3.6 066 */ 067 068public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldElement<T>> 069 extends AdaptiveStepsizeFieldIntegrator<T> 070 implements FieldButcherArrayProvider<T> { 071 072 /** Index of the pre-computed derivative for <i>fsal</i> methods. */ 073 private final int fsal; 074 075 /** Time steps from Butcher array (without the first zero). */ 076 private final T[] c; 077 078 /** Internal weights from Butcher array (without the first empty row). */ 079 private final T[][] a; 080 081 /** External weights for the high order method from Butcher array. */ 082 private final T[] b; 083 084 /** Stepsize control exponent. */ 085 private final T exp; 086 087 /** Safety factor for stepsize control. */ 088 private T safety; 089 090 /** Minimal reduction factor for stepsize control. */ 091 private T minReduction; 092 093 /** Maximal growth factor for stepsize control. */ 094 private T maxGrowth; 095 096 /** Build a Runge-Kutta integrator with the given Butcher array. 097 * @param field field to which the time and state vector elements belong 098 * @param name name of the method 099 * @param fsal index of the pre-computed derivative for <i>fsal</i> methods 100 * or -1 if method is not <i>fsal</i> 101 * @param minStep minimal step (sign is irrelevant, regardless of 102 * integration direction, forward or backward), the last step can 103 * be smaller than this 104 * @param maxStep maximal step (sign is irrelevant, regardless of 105 * integration direction, forward or backward), the last step can 106 * be smaller than this 107 * @param scalAbsoluteTolerance allowed absolute error 108 * @param scalRelativeTolerance allowed relative error 109 */ 110 protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal, 111 final double minStep, final double maxStep, 112 final double scalAbsoluteTolerance, 113 final double scalRelativeTolerance) { 114 115 super(field, name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); 116 117 this.fsal = fsal; 118 this.c = getC(); 119 this.a = getA(); 120 this.b = getB(); 121 122 exp = field.getOne().divide(-getOrder()); 123 124 // set the default values of the algorithm control parameters 125 setSafety(field.getZero().add(0.9)); 126 setMinReduction(field.getZero().add(0.2)); 127 setMaxGrowth(field.getZero().add(10.0)); 128 } 129 130 /** Build a Runge-Kutta integrator with the given Butcher array. 131 * @param field field to which the time and state vector elements belong 132 * @param name name of the method 133 * @param fsal index of the pre-computed derivative for <i>fsal</i> methods 134 * or -1 if method is not <i>fsal</i> 135 * @param minStep minimal step (must be positive even for backward 136 * integration), the last step can be smaller than this 137 * @param maxStep maximal step (must be positive even for backward 138 * integration) 139 * @param vecAbsoluteTolerance allowed absolute error 140 * @param vecRelativeTolerance allowed relative error 141 */ 142 protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal, 143 final double minStep, final double maxStep, 144 final double[] vecAbsoluteTolerance, 145 final double[] vecRelativeTolerance) { 146 147 super(field, name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); 148 149 this.fsal = fsal; 150 this.c = getC(); 151 this.a = getA(); 152 this.b = getB(); 153 154 exp = field.getOne().divide(-getOrder()); 155 156 // set the default values of the algorithm control parameters 157 setSafety(field.getZero().add(0.9)); 158 setMinReduction(field.getZero().add(0.2)); 159 setMaxGrowth(field.getZero().add(10.0)); 160 } 161 162 /** Create a fraction. 163 * @param p numerator 164 * @param q denominator 165 * @return p/q computed in the instance field 166 */ 167 protected T fraction(final int p, final int q) { 168 return getField().getOne().multiply(p).divide(q); 169 } 170 171 /** Create a fraction. 172 * @param p numerator 173 * @param q denominator 174 * @return p/q computed in the instance field 175 */ 176 protected T fraction(final double p, final double q) { 177 return getField().getOne().multiply(p).divide(q); 178 } 179 180 /** Create an interpolator. 181 * @param forward integration direction indicator 182 * @param yDotK slopes at the intermediate points 183 * @param globalPreviousState start of the global step 184 * @param globalCurrentState end of the global step 185 * @param mapper equations mapper for the all equations 186 * @return external weights for the high order method from Butcher array 187 */ 188 protected abstract RungeKuttaFieldStepInterpolator<T> createInterpolator(boolean forward, T[][] yDotK, 189 FieldODEStateAndDerivative<T> globalPreviousState, 190 FieldODEStateAndDerivative<T> globalCurrentState, 191 FieldEquationsMapper<T> mapper); 192 /** Get the order of the method. 193 * @return order of the method 194 */ 195 public abstract int getOrder(); 196 197 /** Get the safety factor for stepsize control. 198 * @return safety factor 199 */ 200 public T getSafety() { 201 return safety; 202 } 203 204 /** Set the safety factor for stepsize control. 205 * @param safety safety factor 206 */ 207 public void setSafety(final T safety) { 208 this.safety = safety; 209 } 210 211 /** {@inheritDoc} */ 212 @Override 213 public FieldODEStateAndDerivative<T> integrate(final FieldExpandableODE<T> equations, 214 final FieldODEState<T> initialState, final T finalTime) 215 throws NumberIsTooSmallException, DimensionMismatchException, 216 MaxCountExceededException, NoBracketingException { 217 218 sanityChecks(initialState, finalTime); 219 final T t0 = initialState.getTime(); 220 final T[] y0 = equations.getMapper().mapState(initialState); 221 setStepStart(initIntegration(equations, t0, y0, finalTime)); 222 final boolean forward = finalTime.subtract(initialState.getTime()).getReal() > 0; 223 224 // create some internal working arrays 225 final int stages = c.length + 1; 226 T[] y = y0; 227 final T[][] yDotK = MathArrays.buildArray(getField(), stages, -1); 228 final T[] yTmp = MathArrays.buildArray(getField(), y0.length); 229 230 // set up integration control objects 231 T hNew = getField().getZero(); 232 boolean firstTime = true; 233 234 // main integration loop 235 setIsLastStep(false); 236 do { 237 238 // iterate over step size, ensuring local normalized error is smaller than 1 239 T error = getField().getZero().add(10); 240 while (error.subtract(1.0).getReal() >= 0) { 241 242 // first stage 243 y = equations.getMapper().mapState(getStepStart()); 244 yDotK[0] = equations.getMapper().mapDerivative(getStepStart()); 245 246 if (firstTime) { 247 final T[] scale = MathArrays.buildArray(getField(), mainSetDimension); 248 if (vecAbsoluteTolerance == null) { 249 for (int i = 0; i < scale.length; ++i) { 250 scale[i] = y[i].abs().multiply(scalRelativeTolerance).add(scalAbsoluteTolerance); 251 } 252 } else { 253 for (int i = 0; i < scale.length; ++i) { 254 scale[i] = y[i].abs().multiply(vecRelativeTolerance[i]).add(vecAbsoluteTolerance[i]); 255 } 256 } 257 hNew = initializeStep(forward, getOrder(), scale, getStepStart(), equations.getMapper()); 258 firstTime = false; 259 } 260 261 setStepSize(hNew); 262 if (forward) { 263 if (getStepStart().getTime().add(getStepSize()).subtract(finalTime).getReal() >= 0) { 264 setStepSize(finalTime.subtract(getStepStart().getTime())); 265 } 266 } else { 267 if (getStepStart().getTime().add(getStepSize()).subtract(finalTime).getReal() <= 0) { 268 setStepSize(finalTime.subtract(getStepStart().getTime())); 269 } 270 } 271 272 // next stages 273 for (int k = 1; k < stages; ++k) { 274 275 for (int j = 0; j < y0.length; ++j) { 276 T sum = yDotK[0][j].multiply(a[k-1][0]); 277 for (int l = 1; l < k; ++l) { 278 sum = sum.add(yDotK[l][j].multiply(a[k-1][l])); 279 } 280 yTmp[j] = y[j].add(getStepSize().multiply(sum)); 281 } 282 283 yDotK[k] = computeDerivatives(getStepStart().getTime().add(getStepSize().multiply(c[k-1])), yTmp); 284 } 285 286 // estimate the state at the end of the step 287 for (int j = 0; j < y0.length; ++j) { 288 T sum = yDotK[0][j].multiply(b[0]); 289 for (int l = 1; l < stages; ++l) { 290 sum = sum.add(yDotK[l][j].multiply(b[l])); 291 } 292 yTmp[j] = y[j].add(getStepSize().multiply(sum)); 293 } 294 295 // estimate the error at the end of the step 296 error = estimateError(yDotK, y, yTmp, getStepSize()); 297 if (error.subtract(1.0).getReal() >= 0) { 298 // reject the step and attempt to reduce error by stepsize control 299 final T factor = RealFieldElement.min(maxGrowth, 300 RealFieldElement.max(minReduction, safety.multiply(error.pow(exp)))); 301 hNew = filterStep(getStepSize().multiply(factor), forward, false); 302 } 303 } 304 final T stepEnd = getStepStart().getTime().add(getStepSize()); 305 final T[] yDotTmp = (fsal >= 0) ? yDotK[fsal] : computeDerivatives(stepEnd, yTmp); 306 final FieldODEStateAndDerivative<T> stateTmp = new FieldODEStateAndDerivative<>(stepEnd, yTmp, yDotTmp); 307 308 // local error is small enough: accept the step, trigger events and step handlers 309 System.arraycopy(yTmp, 0, y, 0, y0.length); 310 setStepStart(acceptStep(createInterpolator(forward, yDotK, getStepStart(), stateTmp, equations.getMapper()), 311 finalTime)); 312 313 if (!isLastStep()) { 314 315 // stepsize control for next step 316 final T factor = RealFieldElement.min(maxGrowth, 317 RealFieldElement.max(minReduction, safety.multiply(error.pow(exp)))); 318 final T scaledH = getStepSize().multiply(factor); 319 final T nextT = getStepStart().getTime().add(scaledH); 320 final boolean nextIsLast = forward ? 321 nextT.subtract(finalTime).getReal() >= 0 : 322 nextT.subtract(finalTime).getReal() <= 0; 323 hNew = filterStep(scaledH, forward, nextIsLast); 324 325 final T filteredNextT = getStepStart().getTime().add(hNew); 326 final boolean filteredNextIsLast = forward ? 327 filteredNextT.subtract(finalTime).getReal() >= 0 : 328 filteredNextT.subtract(finalTime).getReal() <= 0; 329 if (filteredNextIsLast) { 330 hNew = finalTime.subtract(getStepStart().getTime()); 331 } 332 } 333 } while (!isLastStep()); 334 335 final FieldODEStateAndDerivative<T> finalState = getStepStart(); 336 resetInternalState(); 337 return finalState; 338 } 339 340 /** Get the minimal reduction factor for stepsize control. 341 * @return minimal reduction factor 342 */ 343 public T getMinReduction() { 344 return minReduction; 345 } 346 347 /** Set the minimal reduction factor for stepsize control. 348 * @param minReduction minimal reduction factor 349 */ 350 public void setMinReduction(final T minReduction) { 351 this.minReduction = minReduction; 352 } 353 354 /** Get the maximal growth factor for stepsize control. 355 * @return maximal growth factor 356 */ 357 public T getMaxGrowth() { 358 return maxGrowth; 359 } 360 361 /** Set the maximal growth factor for stepsize control. 362 * @param maxGrowth maximal growth factor 363 */ 364 public void setMaxGrowth(final T maxGrowth) { 365 this.maxGrowth = maxGrowth; 366 } 367 368 /** Compute the error ratio. 369 * @param yDotK derivatives computed during the first stages 370 * @param y0 estimate of the step at the start of the step 371 * @param y1 estimate of the step at the end of the step 372 * @param h current step 373 * @return error ratio, greater than 1 if step should be rejected 374 */ 375 protected abstract T estimateError(T[][] yDotK, T[] y0, T[] y1, T h); 376}