You are viewing a plain text version of this content. The canonical link for it is here.
Posted to issues@commons.apache.org by "Andreas mueller (JIRA)" <ji...@apache.org> on 2009/12/25 13:08:29 UTC

[jira] Updated: (MATH-325) Improvement of Romberg extrapolation

     [ https://issues.apache.org/jira/browse/MATH-325?page=com.atlassian.jira.plugin.system.issuetabpanels:all-tabpanel ]

Andreas mueller updated MATH-325:
---------------------------------

    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.

<pre>
	/**
	 * 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);
	}
</pre>


  was:
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);
	}




> Improvement of Romberg extrapolation
> ------------------------------------
>
>                 Key: MATH-325
>                 URL: https://issues.apache.org/jira/browse/MATH-325
>             Project: Commons Math
>          Issue Type: Improvement
>    Affects Versions: 2.0
>            Reporter: Andreas mueller
>             Fix For: 2.1
>
>
> 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.
> <pre>
> 	/**
> 	 * 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);
> 	}
> </pre>

-- 
This message is automatically generated by JIRA.
-
You can reply to this email to add a comment to the issue online.