graggbulirschstoerintegrator.java

Apache的common math数学软件包

            if ((error > 1.0e15) || ((k > 1) && (error > maxError))) {
              // error is too big, we reduce the global step
              hNew   = Math.abs(filterStep(stepSize * stabilityReduction, false));
              reject = true;
              loop   = false;
            } else {

              maxError = Math.max(4 * error, 1.0);

              // compute optimal stepsize for this order
              double exp = 1.0 / (2 * k + 1);
              double fac = stepControl2 / Math.pow(error / stepControl1, exp);
              double pow = Math.pow(stepControl3, exp);
              fac = Math.max(pow / stepControl4, Math.min(1 / pow, fac));
              optimalStep[k]     = Math.abs(filterStep(stepSize * fac, true));
              costPerTimeUnit[k] = costPerStep[k] / optimalStep[k];

              // check convergence
              switch (k - targetIter) {

              case -1 :
                if ((targetIter > 1) && ! previousRejected) {

                  // check if we can stop iterations now
                  if (error <= 1.0) {
                    // convergence have been reached just before targetIter
                    loop = false;
                  } else {
                    // estimate if there is a chance convergence will
                    // be reached on next iteration, using the
                    // asymptotic evolution of error
                    double ratio = ((double) sequence [k] * sequence[k+1]) /
                                   (sequence[0] * sequence[0]);
                    if (error > ratio * ratio) {
                      // we don't expect to converge on next iteration
                      // we reject the step immediately and reduce order
                      reject = true;
                      loop   = false;
                      targetIter = k;
                      if ((targetIter > 1) &&
                          (costPerTimeUnit[targetIter-1] <
                           orderControl1 * costPerTimeUnit[targetIter])) {
                        --targetIter;
                      }
                      hNew = optimalStep[targetIter];
                    }
                  }
                }
                break;

              case 0:
                if (error <= 1.0) {
                  // convergence has been reached exactly at targetIter
                  loop = false;
                } else {
                  // estimate if there is a chance convergence will
                  // be reached on next iteration, using the
                  // asymptotic evolution of error
                  double ratio = ((double) sequence[k+1]) / sequence[0];
                  if (error > ratio * ratio) {
                    // we don't expect to converge on next iteration
                    // we reject the step immediately
                    reject = true;
                    loop = false;
                    if ((targetIter > 1) &&
                        (costPerTimeUnit[targetIter-1] <
                         orderControl1 * costPerTimeUnit[targetIter])) {
                      --targetIter;
                    }
                    hNew = optimalStep[targetIter];
                  }
                }
                break;

              case 1 :
                if (error > 1.0) {
                  reject = true;
                  if ((targetIter > 1) &&
                      (costPerTimeUnit[targetIter-1] <
                       orderControl1 * costPerTimeUnit[targetIter])) {
                    --targetIter;
                  }
                  hNew = optimalStep[targetIter];
                }
                loop = false;
                break;

              default :
                if ((firstTime || lastStep) && (error <= 1.0)) {
                  loop = false;
                }
                break;

              }

            }
          }
        }
      }

      // dense output handling
      double hInt = getMaxStep();
      if (denseOutput && ! reject) {

        // extrapolate state at middle point of the step
        for (int j = 1; j <= k; ++j) {
          extrapolate(0, j, diagonal, yMidDots[0]);
        }

        // derivative at end of step
        equations.computeDerivatives(stepStart + stepSize, y1, yDot1);

        int mu = 2 * k - mudif + 3;

        for (int l = 0; l < mu; ++l) {

          // derivative at middle point of the step
          int l2 = l / 2;
          double factor = Math.pow(0.5 * sequence[l2], l);
          int middleIndex = fk[l2].length / 2;
          for (int i = 0; i < y0.length; ++i) {
            yMidDots[l+1][i] = factor * fk[l2][middleIndex + l][i];
          }
          for (int j = 1; j <= k - l2; ++j) {
            factor = Math.pow(0.5 * sequence[j + l2], l);
            middleIndex = fk[l2+j].length / 2;
            for (int i = 0; i < y0.length; ++i) {
              diagonal[j-1][i] = factor * fk[l2+j][middleIndex+l][i];
            }
            extrapolate(l2, j, diagonal, yMidDots[l+1]);
          }
          for (int i = 0; i < y0.length; ++i) {
            yMidDots[l+1][i] *= stepSize;
          }

          // compute centered differences to evaluate next derivatives
          for (int j = (l + 1) / 2; j <= k; ++j) {
            for (int m = fk[j].length - 1; m >= 2 * (l + 1); --m) {
              for (int i = 0; i < y0.length; ++i) {
                fk[j][m][i] -= fk[j][m-2][i];
              }
            }
          }

        }

        if (mu >= 0) {

          // estimate the dense output coefficients
          GraggBulirschStoerStepInterpolator gbsInterpolator
            = (GraggBulirschStoerStepInterpolator) interpolator;
          gbsInterpolator.computeCoefficients(mu, stepSize);

          if (useInterpolationError) {
            // use the interpolation error to limit stepsize
            double interpError = gbsInterpolator.estimateError(scale);
            hInt = Math.abs(stepSize / Math.max(Math.pow(interpError, 1.0 / (mu+4)),
                                                0.01));
            if (interpError > 10.0) {
              hNew = hInt;
              reject = true;
            }
          }

          // Switching functions handling
          if (!reject) {
            interpolator.storeTime(stepStart + stepSize);
            if (switchesHandler.evaluateStep(interpolator)) {
              reject = true;
              hNew = Math.abs(switchesHandler.getEventTime() - stepStart);
            }
          }

        }

        if (!reject) {
          // we will reuse the slope for the beginning of next step
          firstStepAlreadyComputed = true;
          System.arraycopy(yDot1, 0, yDot0, 0, y0.length);
        }

      }

      if (! reject) {

        // store end of step state
        double nextStep = stepStart + stepSize;
        System.arraycopy(y1, 0, y, 0, y0.length);

        switchesHandler.stepAccepted(nextStep, y);
        if (switchesHandler.stop()) {
          lastStep = true;
        }

        // provide the step data to the step handler
        interpolator.storeTime(nextStep);
        handler.handleStep(interpolator, lastStep);
        stepStart = nextStep;

        if (switchesHandler.reset(stepStart, y) && ! lastStep) {
          // some switching function has triggered changes that
          // invalidate the derivatives, we need to recompute them
          firstStepAlreadyComputed = false;
        }

        int optimalIter;
        if (k == 1) {
          optimalIter = 2;
          if (previousRejected) {
            optimalIter = 1;
          }
        } else if (k <= targetIter) {
          optimalIter = k;
          if (costPerTimeUnit[k-1] < orderControl1 * costPerTimeUnit[k]) {
            optimalIter = k-1;
          } else if (costPerTimeUnit[k] < orderControl2 * costPerTimeUnit[k-1]) {
            optimalIter = Math.min(k+1, sequence.length - 2);
          }
        } else {
          optimalIter = k - 1;
          if ((k > 2) &&
              (costPerTimeUnit[k-2] < orderControl1 * costPerTimeUnit[k-1])) {
            optimalIter = k - 2;
          }
          if (costPerTimeUnit[k] < orderControl2 * costPerTimeUnit[optimalIter]) {
            optimalIter = Math.min(k, sequence.length - 2);
          }
        }

        if (previousRejected) {
          // after a rejected step neither order nor stepsize
          // should increase
          targetIter = Math.min(optimalIter, k);
          hNew = Math.min(Math.abs(stepSize), optimalStep[targetIter]);
        } else {
          // stepsize control
          if (optimalIter <= k) {
            hNew = optimalStep[optimalIter];
          } else {
            if ((k < targetIter) &&
                (costPerTimeUnit[k] < orderControl2 * costPerTimeUnit[k-1])) {
              hNew = filterStep(optimalStep[k] *
                                costPerStep[optimalIter+1] / costPerStep[k],
                                false);
            } else {
              hNew = filterStep(optimalStep[k] *
                                costPerStep[optimalIter] / costPerStep[k],
                                false);
            }
          }

          targetIter = optimalIter;

        }

        newStep = true;

      }

      hNew = Math.min(hNew, hInt);
      if (! forward) {
        hNew = -hNew;
      }

      firstTime = false;

      if (reject) {
        lastStep = false;
        previousRejected = true;
      } else {
        previousRejected = false;
      }

    }

  }

  /** maximal order. */
  private int maxOrder;

  /** step size sequence. */
  private int[] sequence;

  /** overall cost of applying step reduction up to iteration k+1,
   *  in number of calls.
   */
  private int[] costPerStep;

  /** cost per unit step. */
  private double[] costPerTimeUnit;

  /** optimal steps for each order. */
  private double[] optimalStep;

  /** extrapolation coefficients. */
  private double[][] coeff;

  /** stability check enabling parameter. */
  private boolean performTest;

  /** maximal number of checks for each iteration. */
  private int maxChecks;

  /** maximal number of iterations for which checks are performed. */
  private int maxIter;

  /** stepsize reduction factor in case of stability check failure. */
  private double stabilityReduction;

  /** first stepsize control factor. */
  private double stepControl1;

  /** second stepsize control factor. */
  private double stepControl2;

  /** third stepsize control factor. */
  private double stepControl3;

  /** fourth stepsize control factor. */
  private double stepControl4;

  /** first order control factor. */
  private double orderControl1;

  /** second order control factor. */
  private double orderControl2;

  /** dense outpute required. */
  private boolean denseOutput;

  /** use interpolation error in stepsize control. */
  private boolean useInterpolationError;

  /** interpolation order control parameter. */
  private int mudif;

}