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; 19 20 import org.hipparchus.exception.MathIllegalArgumentException; 21 import org.hipparchus.exception.MathIllegalStateException; 22 import org.hipparchus.ode.AbstractIntegrator; 23 import org.hipparchus.ode.ODEState; 24 import org.hipparchus.ode.ODEStateAndDerivative; 25 import org.hipparchus.util.FastMath; 26 27 /** 28 * This abstract class holds the common part of all adaptive 29 * stepsize integrators for Ordinary Differential Equations. 30 * 31 * <p>These algorithms perform integration with stepsize control, which 32 * means the user does not specify the integration step but rather a 33 * tolerance on error. The error threshold is computed as 34 * </p> 35 * <pre> 36 * threshold_i = absTol_i + relTol_i * max (abs (ym), abs (ym+1)) 37 * </pre> 38 * <p> 39 * where absTol_i is the absolute tolerance for component i of the 40 * state vector and relTol_i is the relative tolerance for the same 41 * component. The user can also use only two scalar values absTol and 42 * relTol which will be used for all components. 43 * </p> 44 * <p> 45 * If the Ordinary Differential Equations is an {@link org.hipparchus.ode.ExpandableODE 46 * extended ODE} rather than a {@link 47 * org.hipparchus.ode.OrdinaryDifferentialEquation basic ODE}, then 48 * <em>only</em> the {@link org.hipparchus.ode.ExpandableODE#getPrimary() primary part} 49 * of the state vector is used for stepsize control, not the complete state vector. 50 * </p> 51 * 52 * <p>If the estimated error for ym+1 is such that</p> 53 * <pre> 54 * sqrt((sum (errEst_i / threshold_i)^2 ) / n) < 1 55 * </pre> 56 * 57 * <p>(where n is the main set dimension) then the step is accepted, 58 * otherwise the step is rejected and a new attempt is made with a new 59 * stepsize.</p> 60 * 61 * 62 */ 63 64 public abstract class AdaptiveStepsizeIntegrator 65 extends AbstractIntegrator { 66 67 /** Helper for step size control. */ 68 private StepsizeHelper stepsizeHelper; 69 70 /** Build an integrator with the given stepsize bounds. 71 * The default step handler does nothing. 72 * @param name name of the method 73 * @param minStep minimal step (sign is irrelevant, regardless of 74 * integration direction, forward or backward), the last step can 75 * be smaller than this 76 * @param maxStep maximal step (sign is irrelevant, regardless of 77 * integration direction, forward or backward), the last step can 78 * be smaller than this 79 * @param scalAbsoluteTolerance allowed absolute error 80 * @param scalRelativeTolerance allowed relative error 81 */ 82 public AdaptiveStepsizeIntegrator(final String name, 83 final double minStep, final double maxStep, 84 final double scalAbsoluteTolerance, 85 final double scalRelativeTolerance) { 86 super(name); 87 stepsizeHelper = new StepsizeHelper(minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance); 88 resetInternalState(); 89 } 90 91 /** Build an integrator with the given stepsize bounds. 92 * The default step handler does nothing. 93 * @param name name of the method 94 * @param minStep minimal step (sign is irrelevant, regardless of 95 * integration direction, forward or backward), the last step can 96 * be smaller than this 97 * @param maxStep maximal step (sign is irrelevant, regardless of 98 * integration direction, forward or backward), the last step can 99 * be smaller than this 100 * @param vecAbsoluteTolerance allowed absolute error 101 * @param vecRelativeTolerance allowed relative error 102 */ 103 public AdaptiveStepsizeIntegrator(final String name, 104 final double minStep, final double maxStep, 105 final double[] vecAbsoluteTolerance, 106 final double[] vecRelativeTolerance) { 107 super(name); 108 stepsizeHelper = new StepsizeHelper(minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance); 109 resetInternalState(); 110 } 111 112 /** Set the adaptive step size control parameters. 113 * <p> 114 * A side effect of this method is to also reset the initial 115 * step so it will be automatically computed by the integrator 116 * if {@link #setInitialStepSize(double) setInitialStepSize} 117 * is not called by the user. 118 * </p> 119 * @param minimalStep minimal step (must be positive even for backward 120 * integration), the last step can be smaller than this 121 * @param maximalStep maximal step (must be positive even for backward 122 * integration) 123 * @param absoluteTolerance allowed absolute error 124 * @param relativeTolerance allowed relative error 125 */ 126 public void setStepSizeControl(final double minimalStep, final double maximalStep, 127 final double absoluteTolerance, 128 final double relativeTolerance) { 129 stepsizeHelper = new StepsizeHelper(minimalStep, maximalStep, absoluteTolerance, relativeTolerance); 130 } 131 132 /** Set the adaptive step size control parameters. 133 * <p> 134 * A side effect of this method is to also reset the initial 135 * step so it will be automatically computed by the integrator 136 * if {@link #setInitialStepSize(double) setInitialStepSize} 137 * is not called by the user. 138 * </p> 139 * @param minimalStep minimal step (must be positive even for backward 140 * integration), the last step can be smaller than this 141 * @param maximalStep maximal step (must be positive even for backward 142 * integration) 143 * @param absoluteTolerance allowed absolute error 144 * @param relativeTolerance allowed relative error 145 */ 146 public void setStepSizeControl(final double minimalStep, final double maximalStep, 147 final double[] absoluteTolerance, 148 final double[] relativeTolerance) { 149 stepsizeHelper = new StepsizeHelper(minimalStep, maximalStep, absoluteTolerance, relativeTolerance); 150 } 151 152 /** Get the stepsize helper. 153 * @return stepsize helper 154 * @since 2.0 155 */ 156 protected StepsizeHelper getStepSizeHelper() { 157 return stepsizeHelper; 158 } 159 160 /** Set the initial step size. 161 * <p>This method allows the user to specify an initial positive 162 * step size instead of letting the integrator guess it by 163 * itself. If this method is not called before integration is 164 * started, the initial step size will be estimated by the 165 * integrator.</p> 166 * @param initialStepSize initial step size to use (must be positive even 167 * for backward integration ; providing a negative value or a value 168 * outside of the min/max step interval will lead the integrator to 169 * ignore the value and compute the initial step size by itself) 170 */ 171 public void setInitialStepSize(final double initialStepSize) { 172 stepsizeHelper.setInitialStepSize(initialStepSize); 173 } 174 175 /** {@inheritDoc} */ 176 @Override 177 protected void sanityChecks(final ODEState initialState, final double t) 178 throws MathIllegalArgumentException { 179 super.sanityChecks(initialState, t); 180 stepsizeHelper.setMainSetDimension(initialState.getPrimaryStateDimension()); 181 } 182 183 /** Initialize the integration step. 184 * @param forward forward integration indicator 185 * @param order order of the method 186 * @param scale scaling vector for the state vector (can be shorter than state vector) 187 * @param state0 state at integration start time 188 * @return first integration step 189 * @exception MathIllegalStateException if the number of functions evaluations is exceeded 190 * @exception MathIllegalArgumentException if arrays dimensions do not match equations settings 191 */ 192 public double initializeStep(final boolean forward, final int order, final double[] scale, 193 final ODEStateAndDerivative state0) 194 throws MathIllegalArgumentException, MathIllegalStateException { 195 196 if (stepsizeHelper.getInitialStep() > 0) { 197 // use the user provided value 198 return forward ? stepsizeHelper.getInitialStep() : -stepsizeHelper.getInitialStep(); 199 } 200 201 // very rough first guess : h = 0.01 * ||y/scale|| / ||y'/scale|| 202 // this guess will be used to perform an Euler step 203 final double[] y0 = state0.getCompleteState(); 204 final double[] yDot0 = state0.getCompleteDerivative(); 205 double yOnScale2 = 0; 206 double yDotOnScale2 = 0; 207 for (int j = 0; j < scale.length; ++j) { 208 final double ratio = y0[j] / scale[j]; 209 yOnScale2 += ratio * ratio; 210 final double ratioDot = yDot0[j] / scale[j]; 211 yDotOnScale2 += ratioDot * ratioDot; 212 } 213 214 double h = ((yOnScale2 < 1.0e-10) || (yDotOnScale2 < 1.0e-10)) ? 215 1.0e-6 : (0.01 * FastMath.sqrt(yOnScale2 / yDotOnScale2)); 216 if (h > getMaxStep()) { 217 h = getMaxStep(); 218 } 219 if (! forward) { 220 h = -h; 221 } 222 223 // perform an Euler step using the preceding rough guess 224 final double[] y1 = new double[y0.length]; 225 for (int j = 0; j < y0.length; ++j) { 226 y1[j] = y0[j] + h * yDot0[j]; 227 } 228 final double[] yDot1 = computeDerivatives(state0.getTime() + h, y1); 229 230 // estimate the second derivative of the solution 231 double yDDotOnScale = 0; 232 for (int j = 0; j < scale.length; ++j) { 233 final double ratioDotDot = (yDot1[j] - yDot0[j]) / scale[j]; 234 yDDotOnScale += ratioDotDot * ratioDotDot; 235 } 236 yDDotOnScale = FastMath.sqrt(yDDotOnScale) / h; 237 238 // step size is computed such that 239 // h^order * max (||y'/tol||, ||y''/tol||) = 0.01 240 final double maxInv2 = FastMath.max(FastMath.sqrt(yDotOnScale2), yDDotOnScale); 241 final double h1 = (maxInv2 < 1.0e-15) ? 242 FastMath.max(1.0e-6, 0.001 * FastMath.abs(h)) : 243 FastMath.pow(0.01 / maxInv2, 1.0 / order); 244 h = FastMath.min(100.0 * FastMath.abs(h), h1); 245 h = FastMath.max(h, 1.0e-12 * FastMath.abs(state0.getTime())); // avoids cancellation when computing t1 - t0 246 if (h < getMinStep()) { 247 h = getMinStep(); 248 } 249 if (h > getMaxStep()) { 250 h = getMaxStep(); 251 } 252 if (! forward) { 253 h = -h; 254 } 255 256 return h; 257 258 } 259 260 /** Reset internal state to dummy values. */ 261 protected void resetInternalState() { 262 setStepStart(null); 263 setStepSize(stepsizeHelper.getDummyStepsize()); 264 } 265 266 /** Get the minimal step. 267 * @return minimal step 268 */ 269 public double getMinStep() { 270 return stepsizeHelper.getMinStep(); 271 } 272 273 /** Get the maximal step. 274 * @return maximal step 275 */ 276 public double getMaxStep() { 277 return stepsizeHelper.getMaxStep(); 278 } 279 280 }