You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by er...@apache.org on 2016/05/06 23:29:00 UTC
[math] MATH-1015
Repository: commons-math
Updated Branches:
refs/heads/feature-MATH-1015 [created] 93d35c8f2
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/feature-MATH-1015
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);
+ }
+ }
+}