commons-issues mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Sébastien Brisard (JIRA) <j...@apache.org>
Subject [jira] [Commented] (MATH-797) Single step integrators
Date Tue, 26 Jun 2012 13:45:44 GMT

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13401388#comment-13401388
] 

Sébastien Brisard commented on MATH-797:
----------------------------------------

I would be happy to get involved in this, but certainly do not want to tread on anyone's toes...
I guess it all really depends on how urgently you need this feature (I'm quite busy with testing
the vector implementations, right now, plus summer vacations are approaching...).

My understanding of the above conclusion was that it is better to create a new hierarchy of
classes for Gauss-like integrators. Personally, I'd rather we take some time to discuss the
interface common to these integrators. Implementation could then be quite fast, as I could
easily refactor what I've already written, and provide various implementations (see above),
including tests; also we could refactor what has already been done on your side.

So I would say: OK to get involved (but not tomorrow!), but let's discuss first !

May be to fire up the discussion, I could include the interfaces I had designed a few years
back. It's just to initiate the discussion, though... So here goes. I started by defining
a Gauss integration rule
{code:java}
public interface GaussianRule {
	/**
	 * Copies the abscissae of all integration points in the specified array,
	 * starting at the specified location. This method will generate an
	 * {@link ArrayIndexOutOfBoundsException} if the specified array is not
	 * large enough.
	 * 
	 * @param x
	 *            the array in which the nodes are to be copied
	 * @param pos
	 *            the index of the first cell to be filled
	 */
	void copyNodes(final double[] x, final int pos);

	/**
	 * Copies the weights of all integration points in the specified array,
	 * starting at the specified location. This method will generate an
	 * {@link ArrayIndexOutOfBoundsException} if the specified array is not
	 * large enough.
	 * 
	 * @param w
	 *            the array in which the weights are to be copied
	 * @param pos
	 *            the index of the first cell to be filled
	 */
	void copyWeights(final double[] w, final int pos);

	/**
	 * Returns the abscissa of the <code>i</code>-th integration point.
	 * 
	 * @param i
	 *            the index of the node
	 * @return the abscissa of the node
	 */
	double getNode(final int i);

	/**
	 * Returns the number of integration points of this quadrature rule.
	 * 
	 * @return the number of integration points
	 */
	int getNumNodes();

	/**
	 * Returns the value of the weight of the <code>i</code>-th integration
	 * point.
	 * 
	 * @param i
	 *            the index of the node
	 * @return the weight of the node
	 */
	double getWeight(final int i);

	/**
	 * Returns the weighted integral of the specified function between the
	 * natural bounds of this quadrature rule.
	 * 
	 * @param f
	 *            the univariate real function to be integrated
	 * @return the estimate of the integral
	 */
	double integrate(final UnivariateRealFunction f);

	/**
	 * Returns the maximum degree of polynomials for which integration with this
	 * quadrature rule is exact.
	 * 
	 * @return the maximum degree
	 */
	int getMaximumDegreeForExactness();

	/**
	 * Returns the lower-bound of the integration interval for which this rule
	 * is designed. Appropriate change of the integration variable may allow
	 * integration on a different interval.
	 * 
	 * @return the lower bound of the integration interval
	 */
	double getRangeOfIntegrationLowerBound();

	/**
	 * Returns the upper-bound of the integration interval for which this rule
	 * is designed. Appropriate change of the integration variable may allow
	 * integration on a different interval.
	 * 
	 * @return the upper bound of the integration interval
	 */
	double getRangeOfIntegrationUpperBound();
}
{code}

Then, each integration rule (Gauss-Legendre, Gauss-Chebyshev, etc...) is created with a factory
{code:java}
public interface GaussianRuleFactory {
    /**
     * Returns a new <code>n</code>-point Gaussian quadrature rule. If no rule
     * has been defined for the specified number of integration points, this
     * method returns <code>null</code>.
     * 
     * @param n
     *            the number of integration points
     * @return the quadrature rule
     */
    GaussianRule create(final int n);
}

{code}

Example for Gauss-Legendre ({{@mathml}} tags are taglets I've written to include MathML code
into Javadocs, please do not mind them).
{code:java}
package org.jquadrature.gauss;

import org.changemyname.WeightedPoint1d;

import static org.celestin.constants.DoubleConstants.ONE_HALF;

//@formatter:off
/**
 * <p>
 * This factory creates Gauss-type quadrature rule using Legendre polynomials.
 * Such a quadrature rule allows the exact calculation of integrals
 * {@mathML doc-files/gauss-legendre-01.mml} where {@inlineMathML <mi>f</mi>}
 * is a polynomial of degree lower than or equal to
 * {@mathML doc-files/gauss-legendre-02.mml}, and {@inlineMathML n} is the
 * number of integration points. The lower- and upper-bounds of the natural
 * interval of integration are {@inlineMathML <mo>-</mo><mn>1</mn>}
and 
 * {@inlineMathML <mn>1</mn>} in this implementation.
 * </p>
 * <p>
 * First, the roots of the Legendre polynomial are computed (using bisection)
 * with half-ulp accuracy. Legendre polynomials are evaluated using the
 * recurrence relation (
 * <a href="{@docRoot}/doc-files/references.xhtml#ABRA1964">
 * Abramowitz and Stegun, 1964
 * </a>, 22.7.10) {@mathML doc-files/gauss-legendre-03.mml}
 * </p>
 * <p>
 * Relations 25.4.29 and 22.8.5 from
 * <a href="{@docRoot}/doc-files/references.xhtml#ABRA1964">
 * Abramowitz and Stegun (1964)</a> are then combined to compute the weights
 * {@mathML doc-files/gauss-legendre-04.mml} where {@inlineMathML
 * <msub><mi>x</mi><mi>i</mi></msub>} is the {@inlineMathML
 * <mi>i</mi>}-th node, and {@inlineMathML <msub><mi>w</mi><mi>i</mi></msub>}
 * the corresponding weight. The accuracy of this computation is somewhat worse,
 * that is why all quadrature rules up to order 20 are precomputed using
 * <a href="doc-files/gauss-legendre.mac">this script</a> for the
 * <a href="http://maxima.sourceforge.net"/>Maxima</a> computer algebra system.
 * </p>
 * <p>
 * The root-finding method is based on the property of interlacing of roots of
 * orthogonal polynomials. Therefore, in order to compute the nodes and weights
 * for the {@inlineMathML <mi>n</mi>}-point quadrature rule, the roots of
 * <em>all</em> {@mathML doc-files/gauss-legendre-05.mml} Legendre
 * polynomials must first be computed. These roots are saved for further use.
 * </p>
 * 
 * @author S&eacute;bastien Brisard
 */
//@formatter:on
public class GaussLegendreFactory extends DefaultGaussianRuleFactory {
	/**
	 * The unique instance of this class (singleton pattern).
	 */
	private static GaussianRuleFactory factory;

	/**
	 * The last rule computed.
	 */
	private WeightedPoint1d[] oldRule;

	/**
	 * Creates a new instance of this class.
	 */
	private GaussLegendreFactory() {
		super();

		// 1-point quadrature rule
		oldRule = new WeightedPoint1d[1];
		oldRule[0] = new WeightedPoint1d(0.0E0, 2.0E0);
		addRule(oldRule);

		// 2-point quadrature rule
		oldRule = new WeightedPoint1d[2];
		oldRule[0] = new WeightedPoint1d(-5.7735026918962576451E-1, 1.0E0);
		oldRule[1] = new WeightedPoint1d(5.7735026918962576451E-1, 1.0E0);
		addRule(oldRule);

		// 3-point quadrature rule
		oldRule = new WeightedPoint1d[3];
		oldRule[0] = new WeightedPoint1d(-7.7459666924148337703E-1,
				5.5555555555555555556E-1);
		oldRule[1] = new WeightedPoint1d(0.0E0, 8.8888888888888888889E-1);
		oldRule[2] = new WeightedPoint1d(7.7459666924148337703E-1,
				5.5555555555555555556E-1);
		addRule(oldRule);

		// 4-point quadrature rule
		oldRule = new WeightedPoint1d[4];
		oldRule[0] = new WeightedPoint1d(-8.6113631159405257523E-1,
				3.4785484513745385737E-1);
		oldRule[1] = new WeightedPoint1d(-3.399810435848562648E-1,
				6.5214515486254614263E-1);
		oldRule[2] = new WeightedPoint1d(3.399810435848562648E-1,
				6.5214515486254614263E-1);
		oldRule[3] = new WeightedPoint1d(8.6113631159405257523E-1,
				3.4785484513745385737E-1);
		addRule(oldRule);

		// 5-point quadrature rule
		oldRule = new WeightedPoint1d[5];
		oldRule[0] = new WeightedPoint1d(-9.061798459386639928E-1,
				2.3692688505618908751E-1);
		oldRule[1] = new WeightedPoint1d(-5.3846931010568309104E-1,
				4.7862867049936646804E-1);
		oldRule[2] = new WeightedPoint1d(0.0E0, 5.6888888888888888889E-1);
		oldRule[3] = new WeightedPoint1d(5.3846931010568309104E-1,
				4.7862867049936646804E-1);
		oldRule[4] = new WeightedPoint1d(9.061798459386639928E-1,
				2.3692688505618908751E-1);
		addRule(oldRule);

		// 6-point quadrature rule
		oldRule = new WeightedPoint1d[6];
		oldRule[0] = new WeightedPoint1d(-9.3246951420315202782E-1,
				1.7132449237917034503E-1);
		oldRule[1] = new WeightedPoint1d(-6.6120938646626451366E-1,
				3.6076157304813860757E-1);
		oldRule[2] = new WeightedPoint1d(-2.3861918608319690863E-1,
				4.6791393457269104739E-1);
		oldRule[3] = new WeightedPoint1d(2.3861918608319690863E-1,
				4.6791393457269104739E-1);
		oldRule[4] = new WeightedPoint1d(6.6120938646626451366E-1,
				3.6076157304813860757E-1);
		oldRule[5] = new WeightedPoint1d(9.3246951420315202782E-1,
				1.7132449237917034503E-1);
		addRule(oldRule);

		// 7-point quadrature rule
		oldRule = new WeightedPoint1d[7];
		oldRule[0] = new WeightedPoint1d(-9.4910791234275852452E-1,
				1.2948496616886969328E-1);
		oldRule[1] = new WeightedPoint1d(-7.4153118559939443987E-1,
				2.797053914892766679E-1);
		oldRule[2] = new WeightedPoint1d(-4.0584515137739716691E-1,
				3.8183005050511894495E-1);
		oldRule[3] = new WeightedPoint1d(0.0E0, 4.1795918367346938776E-1);
		oldRule[4] = new WeightedPoint1d(4.0584515137739716691E-1,
				3.8183005050511894495E-1);
		oldRule[5] = new WeightedPoint1d(7.4153118559939443987E-1,
				2.797053914892766679E-1);
		oldRule[6] = new WeightedPoint1d(9.4910791234275852452E-1,
				1.2948496616886969328E-1);
		addRule(oldRule);

		// 8-point quadrature rule
		oldRule = new WeightedPoint1d[8];
		oldRule[0] = new WeightedPoint1d(-9.6028985649753623169E-1,
				1.0122853629037625915E-1);
		oldRule[1] = new WeightedPoint1d(-7.9666647741362673959E-1,
				2.2238103445337447055E-1);
		oldRule[2] = new WeightedPoint1d(-5.2553240991632898582E-1,
				3.1370664587788728734E-1);
		oldRule[3] = new WeightedPoint1d(-1.8343464249564980494E-1,
				3.6268378337836198297E-1);
		oldRule[4] = new WeightedPoint1d(1.8343464249564980494E-1,
				3.6268378337836198297E-1);
		oldRule[5] = new WeightedPoint1d(5.2553240991632898582E-1,
				3.1370664587788728734E-1);
		oldRule[6] = new WeightedPoint1d(7.9666647741362673959E-1,
				2.2238103445337447055E-1);
		oldRule[7] = new WeightedPoint1d(9.6028985649753623169E-1,
				1.0122853629037625915E-1);
		addRule(oldRule);

		// 9-point quadrature rule
		oldRule = new WeightedPoint1d[9];
		oldRule[0] = new WeightedPoint1d(-9.6816023950762608983E-1,
				8.1274388361574411976E-2);
		oldRule[1] = new WeightedPoint1d(-8.360311073266357943E-1,
				1.8064816069485740406E-1);
		oldRule[2] = new WeightedPoint1d(-6.1337143270059039731E-1,
				2.6061069640293546232E-1);
		oldRule[3] = new WeightedPoint1d(-3.2425342340380892904E-1,
				3.1234707704000284007E-1);
		oldRule[4] = new WeightedPoint1d(0.0E0, 3.3023935500125976317E-1);
		oldRule[5] = new WeightedPoint1d(3.2425342340380892904E-1,
				3.1234707704000284007E-1);
		oldRule[6] = new WeightedPoint1d(6.1337143270059039731E-1,
				2.6061069640293546232E-1);
		oldRule[7] = new WeightedPoint1d(8.360311073266357943E-1,
				1.8064816069485740406E-1);
		oldRule[8] = new WeightedPoint1d(9.6816023950762608983E-1,
				8.1274388361574411976E-2);
		addRule(oldRule);

		// 10-point quadrature rule
		oldRule = new WeightedPoint1d[10];
		oldRule[0] = new WeightedPoint1d(-9.7390652851717172008E-1,
				6.6671344308688137597E-2);
		oldRule[1] = new WeightedPoint1d(-8.6506336668898451074E-1,
				1.4945134915058059314E-1);
		oldRule[2] = new WeightedPoint1d(-6.7940956829902440623E-1,
				2.19086362515982044E-1);
		oldRule[3] = new WeightedPoint1d(-4.333953941292471908E-1,
				2.6926671930999635509E-1);
		oldRule[4] = new WeightedPoint1d(-1.4887433898163121088E-1,
				2.9552422471475287018E-1);
		oldRule[5] = new WeightedPoint1d(1.4887433898163121088E-1,
				2.9552422471475287018E-1);
		oldRule[6] = new WeightedPoint1d(4.333953941292471908E-1,
				2.6926671930999635509E-1);
		oldRule[7] = new WeightedPoint1d(6.7940956829902440623E-1,
				2.19086362515982044E-1);
		oldRule[8] = new WeightedPoint1d(8.6506336668898451074E-1,
				1.4945134915058059314E-1);
		oldRule[9] = new WeightedPoint1d(9.7390652851717172008E-1,
				6.6671344308688137597E-2);
		addRule(oldRule);

		// 11-point quadrature rule
		oldRule = new WeightedPoint1d[11];
		oldRule[0] = new WeightedPoint1d(-9.7822865814605699281E-1,
				5.5668567116173666479E-2);
		oldRule[1] = new WeightedPoint1d(-8.8706259976809529907E-1,
				1.2558036946490462464E-1);
		oldRule[2] = new WeightedPoint1d(-7.3015200557404932409E-1,
				1.8629021092773425143E-1);
		oldRule[3] = new WeightedPoint1d(-5.1909612920681181592E-1,
				2.3319376459199047992E-1);
		oldRule[4] = new WeightedPoint1d(-2.6954315595234497233E-1,
				2.6280454451024666218E-1);
		oldRule[5] = new WeightedPoint1d(0.0E0, 2.7292508677790063072E-1);
		oldRule[6] = new WeightedPoint1d(2.6954315595234497233E-1,
				2.6280454451024666218E-1);
		oldRule[7] = new WeightedPoint1d(5.1909612920681181592E-1,
				2.3319376459199047992E-1);
		oldRule[8] = new WeightedPoint1d(7.3015200557404932409E-1,
				1.8629021092773425143E-1);
		oldRule[9] = new WeightedPoint1d(8.8706259976809529907E-1,
				1.2558036946490462464E-1);
		oldRule[10] = new WeightedPoint1d(9.7822865814605699281E-1,
				5.5668567116173666479E-2);
		addRule(oldRule);

		// 12-point quadrature rule
		oldRule = new WeightedPoint1d[12];
		oldRule[0] = new WeightedPoint1d(-9.8156063424671925069E-1,
				4.71753363865118272E-2);
		oldRule[1] = new WeightedPoint1d(-9.0411725637047485668E-1,
				1.0693932599531843096E-1);
		oldRule[2] = new WeightedPoint1d(-7.6990267419430468704E-1,
				1.6007832854334622634E-1);
		oldRule[3] = new WeightedPoint1d(-5.873179542866174473E-1,
				2.0316742672306592175E-1);
		oldRule[4] = new WeightedPoint1d(-3.6783149899818019376E-1,
				2.3349253653835480876E-1);
		oldRule[5] = new WeightedPoint1d(-1.2523340851146891547E-1,
				2.49147045813402785E-1);
		oldRule[6] = new WeightedPoint1d(1.2523340851146891547E-1,
				2.49147045813402785E-1);
		oldRule[7] = new WeightedPoint1d(3.6783149899818019376E-1,
				2.3349253653835480876E-1);
		oldRule[8] = new WeightedPoint1d(5.873179542866174473E-1,
				2.0316742672306592175E-1);
		oldRule[9] = new WeightedPoint1d(7.6990267419430468704E-1,
				1.6007832854334622634E-1);
		oldRule[10] = new WeightedPoint1d(9.0411725637047485668E-1,
				1.0693932599531843096E-1);
		oldRule[11] = new WeightedPoint1d(9.8156063424671925069E-1,
				4.71753363865118272E-2);
		addRule(oldRule);

		// 13-point quadrature rule
		oldRule = new WeightedPoint1d[13];
		oldRule[0] = new WeightedPoint1d(-9.8418305471858814948E-1,
				4.0484004765315879512E-2);
		oldRule[1] = new WeightedPoint1d(-9.1759839922297796521E-1,
				9.2121499837728447916E-2);
		oldRule[2] = new WeightedPoint1d(-8.015780907333099128E-1,
				1.3887351021978723846E-1);
		oldRule[3] = new WeightedPoint1d(-6.4234933944034022064E-1,
				1.7814598076194573828E-1);
		oldRule[4] = new WeightedPoint1d(-4.4849275103644685288E-1,
				2.0781604753688850231E-1);
		oldRule[5] = new WeightedPoint1d(-2.3045831595513479406E-1,
				2.2628318026289723841E-1);
		oldRule[6] = new WeightedPoint1d(0.0E0, 2.3255155323087391019E-1);
		oldRule[7] = new WeightedPoint1d(2.3045831595513479406E-1,
				2.2628318026289723841E-1);
		oldRule[8] = new WeightedPoint1d(4.4849275103644685288E-1,
				2.0781604753688850231E-1);
		oldRule[9] = new WeightedPoint1d(6.4234933944034022064E-1,
				1.7814598076194573828E-1);
		oldRule[10] = new WeightedPoint1d(8.015780907333099128E-1,
				1.3887351021978723846E-1);
		oldRule[11] = new WeightedPoint1d(9.1759839922297796521E-1,
				9.2121499837728447916E-2);
		oldRule[12] = new WeightedPoint1d(9.8418305471858814948E-1,
				4.0484004765315879512E-2);
		addRule(oldRule);

		// 14-point quadrature rule
		oldRule = new WeightedPoint1d[14];
		oldRule[0] = new WeightedPoint1d(-9.8628380869681233884E-1,
				3.5119460331751863032E-2);
		oldRule[1] = new WeightedPoint1d(-9.2843488366357351733E-1,
				8.0158087159760209808E-2);
		oldRule[2] = new WeightedPoint1d(-8.2720131506976499319E-1,
				1.2151857068790318469E-1);
		oldRule[3] = new WeightedPoint1d(-6.8729290481168547015E-1,
				1.5720316715819353457E-1);
		oldRule[4] = new WeightedPoint1d(-5.1524863635815409197E-1,
				1.8553839747793781374E-1);
		oldRule[5] = new WeightedPoint1d(-3.1911236892788976044E-1,
				2.0519846372129560397E-1);
		oldRule[6] = new WeightedPoint1d(-1.0805494870734366207E-1,
				2.152638534631577902E-1);
		oldRule[7] = new WeightedPoint1d(1.0805494870734366207E-1,
				2.152638534631577902E-1);
		oldRule[8] = new WeightedPoint1d(3.1911236892788976044E-1,
				2.0519846372129560397E-1);
		oldRule[9] = new WeightedPoint1d(5.1524863635815409197E-1,
				1.8553839747793781374E-1);
		oldRule[10] = new WeightedPoint1d(6.8729290481168547015E-1,
				1.5720316715819353457E-1);
		oldRule[11] = new WeightedPoint1d(8.2720131506976499319E-1,
				1.2151857068790318469E-1);
		oldRule[12] = new WeightedPoint1d(9.2843488366357351733E-1,
				8.0158087159760209808E-2);
		oldRule[13] = new WeightedPoint1d(9.8628380869681233884E-1,
				3.5119460331751863032E-2);
		addRule(oldRule);

		// 15-point quadrature rule
		oldRule = new WeightedPoint1d[15];
		oldRule[0] = new WeightedPoint1d(-9.8799251802048542849E-1,
				3.0753241996117268357E-2);
		oldRule[1] = new WeightedPoint1d(-9.3727339240070590431E-1,
				7.0366047488108124707E-2);
		oldRule[2] = new WeightedPoint1d(-8.482065834104272162E-1,
				1.0715922046717193501E-1);
		oldRule[3] = new WeightedPoint1d(-7.2441773136017004742E-1,
				1.3957067792615431445E-1);
		oldRule[4] = new WeightedPoint1d(-5.7097217260853884754E-1,
				1.6626920581699393355E-1);
		oldRule[5] = new WeightedPoint1d(-3.941513470775633699E-1,
				1.8616100001556221103E-1);
		oldRule[6] = new WeightedPoint1d(-2.011940939974345223E-1,
				1.9843148532711157646E-1);
		oldRule[7] = new WeightedPoint1d(0.0E0, 2.0257824192556127288E-1);
		oldRule[8] = new WeightedPoint1d(2.011940939974345223E-1,
				1.9843148532711157646E-1);
		oldRule[9] = new WeightedPoint1d(3.941513470775633699E-1,
				1.8616100001556221103E-1);
		oldRule[10] = new WeightedPoint1d(5.7097217260853884754E-1,
				1.6626920581699393355E-1);
		oldRule[11] = new WeightedPoint1d(7.2441773136017004742E-1,
				1.3957067792615431445E-1);
		oldRule[12] = new WeightedPoint1d(8.482065834104272162E-1,
				1.0715922046717193501E-1);
		oldRule[13] = new WeightedPoint1d(9.3727339240070590431E-1,
				7.0366047488108124707E-2);
		oldRule[14] = new WeightedPoint1d(9.8799251802048542849E-1,
				3.0753241996117268357E-2);
		addRule(oldRule);

		// 16-point quadrature rule
		oldRule = new WeightedPoint1d[16];
		oldRule[0] = new WeightedPoint1d(-9.894009349916499326E-1,
				2.7152459411754094847E-2);
		oldRule[1] = new WeightedPoint1d(-9.4457502307323257608E-1,
				6.225352393864789286E-2);
		oldRule[2] = new WeightedPoint1d(-8.6563120238783174388E-1,
				9.5158511682492784808E-2);
		oldRule[3] = new WeightedPoint1d(-7.5540440835500303389E-1,
				1.2462897125553387205E-1);
		oldRule[4] = new WeightedPoint1d(-6.1787624440264374845E-1,
				1.4959598881657673208E-1);
		oldRule[5] = new WeightedPoint1d(-4.5801677765722738634E-1,
				1.6915651939500253819E-1);
		oldRule[6] = new WeightedPoint1d(-2.8160355077925891323E-1,
				1.8260341504492358887E-1);
		oldRule[7] = new WeightedPoint1d(-9.5012509837637440184E-2,
				1.8945061045506849629E-1);
		oldRule[8] = new WeightedPoint1d(9.5012509837637440184E-2,
				1.8945061045506849629E-1);
		oldRule[9] = new WeightedPoint1d(2.8160355077925891323E-1,
				1.8260341504492358887E-1);
		oldRule[10] = new WeightedPoint1d(4.5801677765722738634E-1,
				1.6915651939500253819E-1);
		oldRule[11] = new WeightedPoint1d(6.1787624440264374845E-1,
				1.4959598881657673208E-1);
		oldRule[12] = new WeightedPoint1d(7.5540440835500303389E-1,
				1.2462897125553387205E-1);
		oldRule[13] = new WeightedPoint1d(8.6563120238783174388E-1,
				9.5158511682492784808E-2);
		oldRule[14] = new WeightedPoint1d(9.4457502307323257608E-1,
				6.225352393864789286E-2);
		oldRule[15] = new WeightedPoint1d(9.894009349916499326E-1,
				2.7152459411754094847E-2);
		addRule(oldRule);

		// 17-point quadrature rule
		oldRule = new WeightedPoint1d[17];
		oldRule[0] = new WeightedPoint1d(-9.9057547531441733567E-1,
				2.4148302868547931965E-2);
		oldRule[1] = new WeightedPoint1d(-9.5067552176876776122E-1,
				5.5459529373987201132E-2);
		oldRule[2] = new WeightedPoint1d(-8.8023915372698590213E-1,
				8.5036148317179180882E-2);
		oldRule[3] = new WeightedPoint1d(-7.8151400389680140692E-1,
				1.118838471934039711E-1);
		oldRule[4] = new WeightedPoint1d(-6.5767115921669076585E-1,
				1.3513636846852547329E-1);
		oldRule[5] = new WeightedPoint1d(-5.1269053708647696789E-1,
				1.5404576107681028808E-1);
		oldRule[6] = new WeightedPoint1d(-3.512317634538763153E-1,
				1.6800410215645004451E-1);
		oldRule[7] = new WeightedPoint1d(-1.7848418149584785585E-1,
				1.7656270536699264633E-1);
		oldRule[8] = new WeightedPoint1d(0.0E0, 1.7944647035620652546E-1);
		oldRule[9] = new WeightedPoint1d(1.7848418149584785585E-1,
				1.7656270536699264633E-1);
		oldRule[10] = new WeightedPoint1d(3.512317634538763153E-1,
				1.6800410215645004451E-1);
		oldRule[11] = new WeightedPoint1d(5.1269053708647696789E-1,
				1.5404576107681028808E-1);
		oldRule[12] = new WeightedPoint1d(6.5767115921669076585E-1,
				1.3513636846852547329E-1);
		oldRule[13] = new WeightedPoint1d(7.8151400389680140692E-1,
				1.118838471934039711E-1);
		oldRule[14] = new WeightedPoint1d(8.8023915372698590213E-1,
				8.5036148317179180882E-2);
		oldRule[15] = new WeightedPoint1d(9.5067552176876776122E-1,
				5.5459529373987201132E-2);
		oldRule[16] = new WeightedPoint1d(9.9057547531441733567E-1,
				2.4148302868547931965E-2);
		addRule(oldRule);

		// 18-point quadrature rule
		oldRule = new WeightedPoint1d[18];
		oldRule[0] = new WeightedPoint1d(-9.9156516842093094673E-1,
				2.1616013526483310319E-2);
		oldRule[1] = new WeightedPoint1d(-9.5582394957139775518E-1,
				4.9714548894969796457E-2);
		oldRule[2] = new WeightedPoint1d(-8.926024664975557392E-1,
				7.6425730254889056531E-2);
		oldRule[3] = new WeightedPoint1d(-8.0370495897252311568E-1,
				1.0094204410628716556E-1);
		oldRule[4] = new WeightedPoint1d(-6.9168704306035320787E-1,
				1.2255520671147846019E-1);
		oldRule[5] = new WeightedPoint1d(-5.5977083107394753461E-1,
				1.406429146706506512E-1);
		oldRule[6] = new WeightedPoint1d(-4.1175116146284264604E-1,
				1.5468467512626524493E-1);
		oldRule[7] = new WeightedPoint1d(-2.5188622569150550959E-1,
				1.6427648374583272299E-1);
		oldRule[8] = new WeightedPoint1d(-8.4775013041735301242E-2,
				1.6914238296314359184E-1);
		oldRule[9] = new WeightedPoint1d(8.4775013041735301242E-2,
				1.6914238296314359184E-1);
		oldRule[10] = new WeightedPoint1d(2.5188622569150550959E-1,
				1.6427648374583272299E-1);
		oldRule[11] = new WeightedPoint1d(4.1175116146284264604E-1,
				1.5468467512626524493E-1);
		oldRule[12] = new WeightedPoint1d(5.5977083107394753461E-1,
				1.406429146706506512E-1);
		oldRule[13] = new WeightedPoint1d(6.9168704306035320787E-1,
				1.2255520671147846019E-1);
		oldRule[14] = new WeightedPoint1d(8.0370495897252311568E-1,
				1.0094204410628716556E-1);
		oldRule[15] = new WeightedPoint1d(8.926024664975557392E-1,
				7.6425730254889056531E-2);
		oldRule[16] = new WeightedPoint1d(9.5582394957139775518E-1,
				4.9714548894969796457E-2);
		oldRule[17] = new WeightedPoint1d(9.9156516842093094673E-1,
				2.1616013526483310319E-2);
		addRule(oldRule);

		// 19-point quadrature rule
		oldRule = new WeightedPoint1d[19];
		oldRule[0] = new WeightedPoint1d(-9.9240684384358440319E-1,
				1.9461788229726477028E-2);
		oldRule[1] = new WeightedPoint1d(-9.6020815213483003085E-1,
				4.4814226765699600332E-2);
		oldRule[2] = new WeightedPoint1d(-9.0315590361481790165E-1,
				6.9044542737641226579E-2);
		oldRule[3] = new WeightedPoint1d(-8.2271465653714282498E-1,
				9.1490021622449999465E-2);
		oldRule[4] = new WeightedPoint1d(-7.2096617733522937862E-1,
				1.1156664554733399472E-1);
		oldRule[5] = new WeightedPoint1d(-6.0054530466168102347E-1,
				1.2875396253933622768E-1);
		oldRule[6] = new WeightedPoint1d(-4.6457074137596094572E-1,
				1.4260670217360661178E-1);
		oldRule[7] = new WeightedPoint1d(-3.1656409996362983199E-1,
				1.5276604206585966678E-1);
		oldRule[8] = new WeightedPoint1d(-1.6035864564022537587E-1,
				1.5896884339395434765E-1);
		oldRule[9] = new WeightedPoint1d(0.0E0, 1.6105444984878369598E-1);
		oldRule[10] = new WeightedPoint1d(1.6035864564022537587E-1,
				1.5896884339395434765E-1);
		oldRule[11] = new WeightedPoint1d(3.1656409996362983199E-1,
				1.5276604206585966678E-1);
		oldRule[12] = new WeightedPoint1d(4.6457074137596094572E-1,
				1.4260670217360661178E-1);
		oldRule[13] = new WeightedPoint1d(6.0054530466168102347E-1,
				1.2875396253933622768E-1);
		oldRule[14] = new WeightedPoint1d(7.2096617733522937862E-1,
				1.1156664554733399472E-1);
		oldRule[15] = new WeightedPoint1d(8.2271465653714282498E-1,
				9.1490021622449999465E-2);
		oldRule[16] = new WeightedPoint1d(9.0315590361481790165E-1,
				6.9044542737641226579E-2);
		oldRule[17] = new WeightedPoint1d(9.6020815213483003085E-1,
				4.4814226765699600332E-2);
		oldRule[18] = new WeightedPoint1d(9.9240684384358440319E-1,
				1.9461788229726477028E-2);
		addRule(oldRule);

		// 20-point quadrature rule
		oldRule = new WeightedPoint1d[20];
		oldRule[0] = new WeightedPoint1d(-9.9312859918509492479E-1,
				1.7614007139152118312E-2);
		oldRule[1] = new WeightedPoint1d(-9.6397192727791379127E-1,
				4.0601429800386941329E-2);
		oldRule[2] = new WeightedPoint1d(-9.1223442825132590587E-1,
				6.267204833410906357E-2);
		oldRule[3] = new WeightedPoint1d(-8.391169718222188234E-1,
				8.3276741576704748724E-2);
		oldRule[4] = new WeightedPoint1d(-7.4633190646015079262E-1,
				1.0193011981724043504E-1);
		oldRule[5] = new WeightedPoint1d(-6.3605368072651502545E-1,
				1.1819453196151841731E-1);
		oldRule[6] = new WeightedPoint1d(-5.10867001950827098E-1,
				1.316886384491766269E-1);
		oldRule[7] = new WeightedPoint1d(-3.7370608871541956067E-1,
				1.4209610931838205133E-1);
		oldRule[8] = new WeightedPoint1d(-2.2778585114164507808E-1,
				1.4917298647260374679E-1);
		oldRule[9] = new WeightedPoint1d(-7.6526521133497333755E-2,
				1.527533871307258507E-1);
		oldRule[10] = new WeightedPoint1d(7.6526521133497333755E-2,
				1.527533871307258507E-1);
		oldRule[11] = new WeightedPoint1d(2.2778585114164507808E-1,
				1.4917298647260374679E-1);
		oldRule[12] = new WeightedPoint1d(3.7370608871541956067E-1,
				1.4209610931838205133E-1);
		oldRule[13] = new WeightedPoint1d(5.10867001950827098E-1,
				1.316886384491766269E-1);
		oldRule[14] = new WeightedPoint1d(6.3605368072651502545E-1,
				1.1819453196151841731E-1);
		oldRule[15] = new WeightedPoint1d(7.4633190646015079262E-1,
				1.0193011981724043504E-1);
		oldRule[16] = new WeightedPoint1d(8.391169718222188234E-1,
				8.3276741576704748724E-2);
		oldRule[17] = new WeightedPoint1d(9.1223442825132590587E-1,
				6.267204833410906357E-2);
		oldRule[18] = new WeightedPoint1d(9.6397192727791379127E-1,
				4.0601429800386941329E-2);
		oldRule[19] = new WeightedPoint1d(9.9312859918509492479E-1,
				1.7614007139152118312E-2);
		addRule(oldRule);
	}

	/**
	 * Returns an instance of this class.
	 * 
	 * @return an instance
	 */
	public static GaussianRuleFactory getInstance() {
		if (factory == null) {
			factory = new GaussLegendreFactory();
		}
		return factory;
	}

	@Override
	public GaussianRule create(final int n) {
		while (oldRule.length < n) {
			incrementOrder();
		}
		return super.create(n);
	}

	/**
	 * Computes the nodes and weights of the next Gaussian quadrature rule.
	 */
	private void incrementOrder() {
		// Degree of the last polynomial computed
		final int n = oldRule.length;
		final WeightedPoint1d[] newRule = new WeightedPoint1d[n + 1];
		// Lower-bound of the interval, and corresponding values of P[j-1],
		// P[j], P[j+1].
		double a, pma, pa, ppa;
		// Upper-bound of the interval, and corresponding values of P[j-1],
		// P[j], P[j+1].
		double b, pmb, pb, ppb;
		// Middle of the interval, and corresponding values of P[j-1], P[j],
		// and P[j+1].
		double x, pm, p, pp;
		final int imax = (n + 1) / 2;
		// Find i-th root of P[n+1] by bracketing.
		for (int i = 0; i < imax; i++) {
			if (i == 0) {
				a = -1.;
			} else {
				a = oldRule[i - 1].getX();
			}
			if (i == n) {
				b = 1.;
			} else {
				b = oldRule[i].getX();
			}
			pma = 1.;
			pa = a;
			pmb = 1.;
			pb = b;
			for (int j = 1; j <= n; j++) {
				// Compute P[j+1](a) and P[j+1](b)
				ppa = ((2 * j + 1) * a * pa - j * pma) / (j + 1);
				ppb = ((2 * j + 1) * b * pb - j * pmb) / (j + 1);
				pma = pa;
				pa = ppa;
				pmb = pb;
				pb = ppb;
			}
			// Now pa = P[n+1](a), and pma = P[n](a) (same holds for b).
			x = ONE_HALF * (a + b);
			pm = 1.;
			p = x;
			boolean done = false;
			while (!done) {
				done = b - a <= Math.ulp(x);
				pm = 1.;
				p = x;
				for (int j = 1; j <= n; j++) {
					// Compute P[j+1](c)
					pp = ((2 * j + 1) * x * p - j * pm) / (j + 1);
					pm = p;
					p = pp;
				}
				// Now pc = P[n+1](c) and pmc = P[n](c).
				if (!done) {
					if (pa * p < 0.) {
						b = x;
						pmb = pm;
						pb = p;
					} else {
						a = x;
						pma = pm;
						pa = p;
					}
					x = ONE_HALF * (a + b);
				}
			}
			double dummy = (n + 1) * (pm - x * p);
			double w = 2. * (1. - x * x) / (dummy * dummy);
			newRule[i] = new WeightedPoint1d(x, w);
			newRule[n - i] = new WeightedPoint1d(-x, w);
		}
		// If (n+1) is odd, 0 is a root.
		if (n % 2 == 0) {
			pm = 1.;
			for (int j = 1; j < n; j += 2) {
				pm = -j * pm / (j + 1.);
			}
			double dummy = (n + 1) * pm;
			double w = 2. / (dummy * dummy);
			newRule[imax] = new WeightedPoint1d(0., w);
		}
		addRule(newRule);
		oldRule = newRule;
	}
}
{code}

It should be noted that Gauss-Legendre rules of any order can be computed and stored. Also,
I chose the midpoint rule instead of Newton in order to guarantee the accuracy of the result.

I would be happy to read your thoughts on all that.
Sébastien
                
> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so
the large interval has to be subdivided into many subintervals. CM does the partition, and
performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval
is required. Those use-cases could benefit from the efficiency gain of not performing a convergence
check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator"
would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira

       

Mime
View raw message