MATH-1015
Gauss-Laguerre quadrature scheme. Thanks to Thomas Neidhart.
This commit is contained in:
parent
f695c9ce35
commit
b77184a218
src
main/java/org/apache/commons/math4/analysis/integration/gauss
test/java/org/apache/commons/math4/analysis/integration/gauss
|
@ -35,6 +35,27 @@ public class GaussIntegratorFactory {
|
||||||
private final BaseRuleFactory<BigDecimal> legendreHighPrecision = new LegendreHighPrecisionRuleFactory();
|
private final BaseRuleFactory<BigDecimal> legendreHighPrecision = new LegendreHighPrecisionRuleFactory();
|
||||||
/** Generator of Gauss-Hermite integrators. */
|
/** Generator of Gauss-Hermite integrators. */
|
||||||
private final BaseRuleFactory<Double> hermite = new HermiteRuleFactory();
|
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.
|
* Creates a Gauss-Legendre integrator of the given order.
|
||||||
|
|
|
@ -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;
|
||||||
|
}
|
||||||
|
}
|
|
@ -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);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
Loading…
Reference in New Issue