Uploaded image for project: 'Commons Math'
  1. Commons Math
  2. MATH-325

Improvement of Romberg extrapolation

    XMLWordPrintableJSON

Details

    • Improvement
    • Status: Closed
    • Major
    • Resolution: Fixed
    • 2.0
    • 2.1
    • None

    Description

      One can use a one-dimensional array (instead of Romberg's tableau) for extrapolating subsequent values.
      Please have a look at following code fragments (which I've taken form the class RombergExtrapolator of
      my MathLibrary). Feel free to use this code.

      /**

      • Default number of maximal extrapolation steps.
        */
        public static int DEF_MAXIMAL_EXTRAPOLATION_COUNT = 8;

      /**

      • The approximation order. <br>
      • Assume that f(h) is approximated by a function a(h), so that f(h) = a(h) +
      • O(h<sup>p</sup>). We say that p is the approximation order.
        */
        private int approximationOrder;
        private int extrapolationCount = 0;
        private double prevResult;

      /**

      • The estimate and tolerance may be used to deside wether to finalize the
      • iteration process (|estimate| < tolerance).
        */

      /** Holds the current estimated error. */
      private double estimate;
      /** Holds the current reached tolerance. */
      private double tolerance;

      private double result[] = new double[DEF_MAXIMAL_EXTRAPOLATION_COUNT + 1];;

      /**

      • Set the maximal number of subsequent extrapolation steps.
      • @param maximalExtrapolationCount
      • maximal extrapolation steps
        */
        public void setMaximalExtrapolationCount(int maximalExtrapolationCount) { result = new double[maximalExtrapolationCount + 1]; }

      /**

      • Extrapolate a sequence of values by means of Romberg's algorithm.
      • Therefore a polynomial of degree maximalExtraploationCount
      • is used. Calculates the current estimate and tolerance using the
      • approximation order.
      • @param value
      • value to extrapolate
      • @return extrapolated value
        */
        public double extrapolate(double value)
        {
        if (extrapolationCount == 0) { // first estimate estimate = value; tolerance = -1.0; prevResult = 0; }

      int i, m, m1 = idx(extrapolationCount);
      long k = (1 << approximationOrder);
      int imin = Math.max(0, extrapolationCount - (result.length - 1));

      result[m1] = value;

      for (i = extrapolationCount - 1; i >= imin; i--)

      { m = idx(i); m1 = idx(i + 1); result[m] = (k * result[m1] - result[m]) / (k - 1); k <<= approximationOrder; }

      m1 = idx(i + 1);
      estimate = result[m1] - prevResult;
      tolerance = Math.abs(result[m1]) * relativeAccuracy + absoluteAccuracy;

      prevResult = result[m1];

      extrapolationCount++;

      return result[m1];
      }

      /**

      • Ring buffer index
        */
        private int idx(int i) { return (i % result.length); }

      Attachments

        Activity

          People

            Unassigned Unassigned
            alpha Andreas mueller
            Votes:
            0 Vote for this issue
            Watchers:
            0 Start watching this issue

            Dates

              Created:
              Updated:
              Resolved:

              Time Tracking

                Estimated:
                Original Estimate - 2h
                2h
                Remaining:
                Remaining Estimate - 2h
                2h
                Logged:
                Time Spent - Not Specified
                Not Specified