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
18package org.apache.commons.math.ode.nonstiff;
19
20import org.apache.commons.math.exception.util.LocalizedFormats;
21import org.apache.commons.math.ode.AbstractIntegrator;
22import org.apache.commons.math.ode.DerivativeException;
23import org.apache.commons.math.ode.ExtendedFirstOrderDifferentialEquations;
24import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
25import org.apache.commons.math.ode.IntegratorException;
26import org.apache.commons.math.util.FastMath;
27
28/**
29 * This abstract class holds the common part of all adaptive
30 * stepsize integrators for Ordinary Differential Equations.
31 *
32 * <p>These algorithms perform integration with stepsize control, which
33 * means the user does not specify the integration step but rather a
34 * tolerance on error. The error threshold is computed as
35 * <pre>
36 * threshold_i = absTol_i + relTol_i * max (abs (ym), abs (ym+1))
37 * </pre>
38 * where absTol_i is the absolute tolerance for component i of the
39 * state vector and relTol_i is the relative tolerance for the same
40 * component. The user can also use only two scalar values absTol and
41 * relTol which will be used for all components.
42 * </p>
43 *
44 * <p>If the Ordinary Differential Equations is an {@link ExtendedFirstOrderDifferentialEquations
45 * extended ODE} rather than a {@link FirstOrderDifferentialEquations basic ODE},
46 * then <em>only</em> the {@link ExtendedFirstOrderDifferentialEquations#getMainSetDimension()
47 * main set} part of the state vector is used for stepsize control, not the complete
48 * state vector.
49 * </p>
50 *
51 * <p>If the estimated error for ym+1 is such that
52 * <pre>
53 * sqrt((sum (errEst_i / threshold_i)^2 ) / n) < 1
54 * </pre>
55 *
56 * (where n is the main set dimension) then the step is accepted,
57 * otherwise the step is rejected and a new attempt is made with a new
58 * stepsize.</p>
59 *
60 * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $
61 * @since 1.2
62 *
63 */
64
65public abstract class AdaptiveStepsizeIntegrator
66  extends AbstractIntegrator {
67
68    /** Allowed absolute scalar error. */
69    protected final double scalAbsoluteTolerance;
70
71    /** Allowed relative scalar error. */
72    protected final double scalRelativeTolerance;
73
74    /** Allowed absolute vectorial error. */
75    protected final double[] vecAbsoluteTolerance;
76
77    /** Allowed relative vectorial error. */
78    protected final double[] vecRelativeTolerance;
79
80    /** Main set dimension. */
81    protected int mainSetDimension;
82
83    /** User supplied initial step. */
84    private double initialStep;
85
86    /** Minimal step. */
87    private final double minStep;
88
89    /** Maximal step. */
90    private final double maxStep;
91
92  /** Build an integrator with the given stepsize bounds.
93   * The default step handler does nothing.
94   * @param name name of the method
95   * @param minStep minimal step (must be positive even for backward
96   * integration), the last step can be smaller than this
97   * @param maxStep maximal step (must be positive even for backward
98   * integration)
99   * @param scalAbsoluteTolerance allowed absolute error
100   * @param scalRelativeTolerance allowed relative error
101   */
102  public AdaptiveStepsizeIntegrator(final String name,
103                                    final double minStep, final double maxStep,
104                                    final double scalAbsoluteTolerance,
105                                    final double scalRelativeTolerance) {
106
107    super(name);
108
109    this.minStep     = FastMath.abs(minStep);
110    this.maxStep     = FastMath.abs(maxStep);
111    this.initialStep = -1.0;
112
113    this.scalAbsoluteTolerance = scalAbsoluteTolerance;
114    this.scalRelativeTolerance = scalRelativeTolerance;
115    this.vecAbsoluteTolerance  = null;
116    this.vecRelativeTolerance  = null;
117
118    resetInternalState();
119
120  }
121
122  /** Build an integrator with the given stepsize bounds.
123   * The default step handler does nothing.
124   * @param name name of the method
125   * @param minStep minimal step (must be positive even for backward
126   * integration), the last step can be smaller than this
127   * @param maxStep maximal step (must be positive even for backward
128   * integration)
129   * @param vecAbsoluteTolerance allowed absolute error
130   * @param vecRelativeTolerance allowed relative error
131   */
132  public AdaptiveStepsizeIntegrator(final String name,
133                                    final double minStep, final double maxStep,
134                                    final double[] vecAbsoluteTolerance,
135                                    final double[] vecRelativeTolerance) {
136
137    super(name);
138
139    this.minStep     = minStep;
140    this.maxStep     = maxStep;
141    this.initialStep = -1.0;
142
143    this.scalAbsoluteTolerance = 0;
144    this.scalRelativeTolerance = 0;
145    this.vecAbsoluteTolerance  = vecAbsoluteTolerance.clone();
146    this.vecRelativeTolerance  = vecRelativeTolerance.clone();
147
148    resetInternalState();
149
150  }
151
152  /** Set the initial step size.
153   * <p>This method allows the user to specify an initial positive
154   * step size instead of letting the integrator guess it by
155   * itself. If this method is not called before integration is
156   * started, the initial step size will be estimated by the
157   * integrator.</p>
158   * @param initialStepSize initial step size to use (must be positive even
159   * for backward integration ; providing a negative value or a value
160   * outside of the min/max step interval will lead the integrator to
161   * ignore the value and compute the initial step size by itself)
162   */
163  public void setInitialStepSize(final double initialStepSize) {
164    if ((initialStepSize < minStep) || (initialStepSize > maxStep)) {
165      initialStep = -1.0;
166    } else {
167      initialStep = initialStepSize;
168    }
169  }
170
171  /** Perform some sanity checks on the integration parameters.
172   * @param equations differential equations set
173   * @param t0 start time
174   * @param y0 state vector at t0
175   * @param t target time for the integration
176   * @param y placeholder where to put the state vector
177   * @exception IntegratorException if some inconsistency is detected
178   */
179  @Override
180  protected void sanityChecks(final FirstOrderDifferentialEquations equations,
181                              final double t0, final double[] y0,
182                              final double t, final double[] y)
183      throws IntegratorException {
184
185      super.sanityChecks(equations, t0, y0, t, y);
186
187      if (equations instanceof ExtendedFirstOrderDifferentialEquations) {
188          mainSetDimension = ((ExtendedFirstOrderDifferentialEquations) equations).getMainSetDimension();
189      } else {
190          mainSetDimension = equations.getDimension();
191      }
192
193      if ((vecAbsoluteTolerance != null) && (vecAbsoluteTolerance.length != mainSetDimension)) {
194          throw new IntegratorException(
195                  LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, mainSetDimension, vecAbsoluteTolerance.length);
196      }
197
198      if ((vecRelativeTolerance != null) && (vecRelativeTolerance.length != mainSetDimension)) {
199          throw new IntegratorException(
200                  LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, mainSetDimension, vecRelativeTolerance.length);
201      }
202
203  }
204
205  /** Initialize the integration step.
206   * @param equations differential equations set
207   * @param forward forward integration indicator
208   * @param order order of the method
209   * @param scale scaling vector for the state vector (can be shorter than state vector)
210   * @param t0 start time
211   * @param y0 state vector at t0
212   * @param yDot0 first time derivative of y0
213   * @param y1 work array for a state vector
214   * @param yDot1 work array for the first time derivative of y1
215   * @return first integration step
216   * @exception DerivativeException this exception is propagated to
217   * the caller if the underlying user function triggers one
218   */
219  public double initializeStep(final FirstOrderDifferentialEquations equations,
220                               final boolean forward, final int order, final double[] scale,
221                               final double t0, final double[] y0, final double[] yDot0,
222                               final double[] y1, final double[] yDot1)
223      throws DerivativeException {
224
225    if (initialStep > 0) {
226      // use the user provided value
227      return forward ? initialStep : -initialStep;
228    }
229
230    // very rough first guess : h = 0.01 * ||y/scale|| / ||y'/scale||
231    // this guess will be used to perform an Euler step
232    double ratio;
233    double yOnScale2 = 0;
234    double yDotOnScale2 = 0;
235    for (int j = 0; j < scale.length; ++j) {
236      ratio         = y0[j] / scale[j];
237      yOnScale2    += ratio * ratio;
238      ratio         = yDot0[j] / scale[j];
239      yDotOnScale2 += ratio * ratio;
240    }
241
242    double h = ((yOnScale2 < 1.0e-10) || (yDotOnScale2 < 1.0e-10)) ?
243               1.0e-6 : (0.01 * FastMath.sqrt(yOnScale2 / yDotOnScale2));
244    if (! forward) {
245      h = -h;
246    }
247
248    // perform an Euler step using the preceding rough guess
249    for (int j = 0; j < y0.length; ++j) {
250      y1[j] = y0[j] + h * yDot0[j];
251    }
252    computeDerivatives(t0 + h, y1, yDot1);
253
254    // estimate the second derivative of the solution
255    double yDDotOnScale = 0;
256    for (int j = 0; j < scale.length; ++j) {
257      ratio         = (yDot1[j] - yDot0[j]) / scale[j];
258      yDDotOnScale += ratio * ratio;
259    }
260    yDDotOnScale = FastMath.sqrt(yDDotOnScale) / h;
261
262    // step size is computed such that
263    // h^order * max (||y'/tol||, ||y''/tol||) = 0.01
264    final double maxInv2 = FastMath.max(FastMath.sqrt(yDotOnScale2), yDDotOnScale);
265    final double h1 = (maxInv2 < 1.0e-15) ?
266                      FastMath.max(1.0e-6, 0.001 * FastMath.abs(h)) :
267                      FastMath.pow(0.01 / maxInv2, 1.0 / order);
268    h = FastMath.min(100.0 * FastMath.abs(h), h1);
269    h = FastMath.max(h, 1.0e-12 * FastMath.abs(t0));  // avoids cancellation when computing t1 - t0
270    if (h < getMinStep()) {
271      h = getMinStep();
272    }
273    if (h > getMaxStep()) {
274      h = getMaxStep();
275    }
276    if (! forward) {
277      h = -h;
278    }
279
280    return h;
281
282  }
283
284  /** Filter the integration step.
285   * @param h signed step
286   * @param forward forward integration indicator
287   * @param acceptSmall if true, steps smaller than the minimal value
288   * are silently increased up to this value, if false such small
289   * steps generate an exception
290   * @return a bounded integration step (h if no bound is reach, or a bounded value)
291   * @exception IntegratorException if the step is too small and acceptSmall is false
292   */
293  protected double filterStep(final double h, final boolean forward, final boolean acceptSmall)
294    throws IntegratorException {
295
296      double filteredH = h;
297      if (FastMath.abs(h) < minStep) {
298          if (acceptSmall) {
299              filteredH = forward ? minStep : -minStep;
300          } else {
301              throw new IntegratorException(
302                      LocalizedFormats.MINIMAL_STEPSIZE_REACHED_DURING_INTEGRATION,
303                      minStep, FastMath.abs(h));
304          }
305      }
306
307      if (filteredH > maxStep) {
308          filteredH = maxStep;
309      } else if (filteredH < -maxStep) {
310          filteredH = -maxStep;
311      }
312
313      return filteredH;
314
315  }
316
317  /** {@inheritDoc} */
318  public abstract double integrate (FirstOrderDifferentialEquations equations,
319                                    double t0, double[] y0,
320                                    double t, double[] y)
321    throws DerivativeException, IntegratorException;
322
323  /** {@inheritDoc} */
324  @Override
325  public double getCurrentStepStart() {
326    return stepStart;
327  }
328
329  /** Reset internal state to dummy values. */
330  protected void resetInternalState() {
331    stepStart = Double.NaN;
332    stepSize  = FastMath.sqrt(minStep * maxStep);
333  }
334
335  /** Get the minimal step.
336   * @return minimal step
337   */
338  public double getMinStep() {
339    return minStep;
340  }
341
342  /** Get the maximal step.
343   * @return maximal step
344   */
345  public double getMaxStep() {
346    return maxStep;
347  }
348
349}
350