MATH-1015 Gauss-Laguerre quadrature scheme. Thanks to Thomas Neidhart.
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/93d35c8f Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/93d35c8f Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/93d35c8f Branch: refs/heads/develop Commit: 93d35c8f291b8953f8460c0eb02a51cc6d6cac22 Parents: cbae75b Author: Gilles <er...@apache.org> Authored: Sat May 7 01:25:42 2016 +0200 Committer: Gilles <er...@apache.org> Committed: Sat May 7 01:25:42 2016 +0200 ---------------------------------------------------------------------- .../gauss/GaussIntegratorFactory.java | 21 +++++ .../integration/gauss/LaguerreRuleFactory.java | 84 ++++++++++++++++++++ .../integration/gauss/LaguerreTest.java | 50 ++++++++++++ 3 files changed, 155 insertions(+) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-math/blob/93d35c8f/src/main/java/org/apache/commons/math4/analysis/integration/gauss/GaussIntegratorFactory.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/analysis/integration/gauss/GaussIntegratorFactory.java b/src/main/java/org/apache/commons/math4/analysis/integration/gauss/GaussIntegratorFactory.java index d15d6ab..4dacc67 100644 --- a/src/main/java/org/apache/commons/math4/analysis/integration/gauss/GaussIntegratorFactory.java +++ b/src/main/java/org/apache/commons/math4/analysis/integration/gauss/GaussIntegratorFactory.java @@ -35,6 +35,27 @@ public class GaussIntegratorFactory { private final BaseRuleFactory<BigDecimal> legendreHighPrecision = new LegendreHighPrecisionRuleFactory(); /** Generator of Gauss-Hermite integrators. */ private final BaseRuleFactory<Double> hermite = new HermiteRuleFactory(); + /** Generator of Gauss-Laguerre integrators. */ + private final BaseRuleFactory<Double> laguerre = new LaguerreRuleFactory(); + + /** + * Creates a Gauss-Laguerre integrator of the given order. + * The call to the + * {@link GaussIntegrator#integrate(org.apache.commons.math4.analysis.UnivariateFunction) + * integrate} method will perform an integration on the interval + * \([0, +\infty)\): the computed value is the improper integral of + * \(e^{-x} f(x)\) + * where \(f(x)\) is the function passed to the + * {@link SymmetricGaussIntegrator#integrate(org.apache.commons.math4.analysis.UnivariateFunction) + * integrate} method. + * + * @param numberOfPoints Order of the integration rule. + * @return a Gauss-Legendre integrator. + * @since 4.0 + */ + public GaussIntegrator laguerre(int numberOfPoints) { + return new GaussIntegrator(getRule(laguerre, numberOfPoints)); + } /** * Creates a Gauss-Legendre integrator of the given order. http://git-wip-us.apache.org/repos/asf/commons-math/blob/93d35c8f/src/main/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreRuleFactory.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreRuleFactory.java b/src/main/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreRuleFactory.java new file mode 100644 index 0000000..220c460 --- /dev/null +++ b/src/main/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreRuleFactory.java @@ -0,0 +1,84 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math4.analysis.integration.gauss; + +import java.util.Arrays; + +import org.apache.commons.math4.analysis.polynomials.PolynomialFunction; +import org.apache.commons.math4.analysis.polynomials.PolynomialsUtils; +import org.apache.commons.math4.exception.DimensionMismatchException; +import org.apache.commons.math4.linear.EigenDecomposition; +import org.apache.commons.math4.linear.MatrixUtils; +import org.apache.commons.math4.linear.RealMatrix; +import org.apache.commons.math4.util.Pair; + +/** + * Factory that creates Gauss-type quadrature rule using Laguerre polynomials. + * + * @see <a href="http://en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature";>Gauss-Laguerre quadrature (Wikipedia)</a> + * @since 4.0 + */ +public class LaguerreRuleFactory extends BaseRuleFactory<Double> { + /** {@inheritDoc} */ + @Override + protected Pair<Double[], Double[]> computeRule(int numberOfPoints) + throws DimensionMismatchException { + + final RealMatrix companionMatrix = companionMatrix(numberOfPoints); + final EigenDecomposition eigen = new EigenDecomposition(companionMatrix); + final double[] roots = eigen.getRealEigenvalues(); + Arrays.sort(roots); + + final Double[] points = new Double[numberOfPoints]; + final Double[] weights = new Double[numberOfPoints]; + + final int n1 = numberOfPoints + 1; + final long n1Squared = n1 * (long) n1; + final PolynomialFunction laguerreN1 = PolynomialsUtils.createLaguerrePolynomial(n1); + for (int i = 0; i < numberOfPoints; i++) { + final double xi = roots[i]; + points[i] = xi; + + final double val = laguerreN1.value(xi); + weights[i] = xi / n1Squared / (val * val); + } + + return new Pair<Double[], Double[]>(points, weights); + } + + /** + * @param degree Matrix dimension. + * @return a square matrix. + */ + private RealMatrix companionMatrix(final int degree) { + final RealMatrix c = MatrixUtils.createRealMatrix(degree, degree); + + for (int i = 0; i < degree; i++) { + c.setEntry(i, i, 2 * i + 1); + if (i + 1 < degree) { + // subdiagonal + c.setEntry(i+1, i, -(i + 1)); + } + if (i - 1 >= 0) { + // superdiagonal + c.setEntry(i-1, i, -i); + } + } + + return c; + } +} http://git-wip-us.apache.org/repos/asf/commons-math/blob/93d35c8f/src/test/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreTest.java ---------------------------------------------------------------------- diff --git a/src/test/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreTest.java b/src/test/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreTest.java new file mode 100644 index 0000000..87b502b --- /dev/null +++ b/src/test/java/org/apache/commons/math4/analysis/integration/gauss/LaguerreTest.java @@ -0,0 +1,50 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math4.analysis.integration.gauss; + +import org.apache.commons.math4.analysis.UnivariateFunction; +import org.apache.commons.math4.special.Gamma; +import org.apache.commons.math4.util.FastMath; +import org.junit.Assert; +import org.junit.Test; + +/** + * Test of the {@link LaguerreRuleFactory}. + */ +public class LaguerreTest { + private static final GaussIntegratorFactory factory = new GaussIntegratorFactory(); + + @Test + public void testGamma() { + final double tol = 1e-13; + + for (int i = 2; i < 10; i += 1) { + final double t = i; + + final UnivariateFunction f = new UnivariateFunction() { + @Override + public double value(double x) { + return FastMath.pow(x, t - 1); + } + }; + + final GaussIntegrator integrator = factory.laguerre(7); + final double s = integrator.integrate(f); + Assert.assertEquals(1d, Gamma.gamma(t) / s, tol); + } + } +}
