diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math/analysis/integration/LegendreGaussIntegrator.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math/analysis/integration/LegendreGaussIntegrator.java | 236 |
1 files changed, 236 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math/analysis/integration/LegendreGaussIntegrator.java b/src/main/java/org/apache/commons/math/analysis/integration/LegendreGaussIntegrator.java new file mode 100644 index 0000000..db6b76c --- /dev/null +++ b/src/main/java/org/apache/commons/math/analysis/integration/LegendreGaussIntegrator.java @@ -0,0 +1,236 @@ +/* + * 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.math.analysis.integration; + +import org.apache.commons.math.ConvergenceException; +import org.apache.commons.math.FunctionEvaluationException; +import org.apache.commons.math.MathRuntimeException; +import org.apache.commons.math.MaxIterationsExceededException; +import org.apache.commons.math.analysis.UnivariateRealFunction; +import org.apache.commons.math.exception.util.LocalizedFormats; +import org.apache.commons.math.util.FastMath; + +/** + * Implements the <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html"> + * Legendre-Gauss</a> quadrature formula. + * <p> + * Legendre-Gauss integrators are efficient integrators that can + * accurately integrate functions with few functions evaluations. A + * Legendre-Gauss integrator using an n-points quadrature formula can + * integrate exactly 2n-1 degree polynomials. + * </p> + * <p> + * These integrators evaluate the function on n carefully chosen + * abscissas in each step interval (mapped to the canonical [-1 1] interval). + * The evaluation abscissas are not evenly spaced and none of them are + * at the interval endpoints. This implies the function integrated can be + * undefined at integration interval endpoints. + * </p> + * <p> + * The evaluation abscissas x<sub>i</sub> are the roots of the degree n + * Legendre polynomial. The weights a<sub>i</sub> of the quadrature formula + * integrals from -1 to +1 ∫ Li<sup>2</sup> where Li (x) = + * ∏ (x-x<sub>k</sub>)/(x<sub>i</sub>-x<sub>k</sub>) for k != i. + * </p> + * <p> + * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 févr. 2011) $ + * @since 1.2 + */ + +public class LegendreGaussIntegrator extends UnivariateRealIntegratorImpl { + + /** Abscissas for the 2 points method. */ + private static final double[] ABSCISSAS_2 = { + -1.0 / FastMath.sqrt(3.0), + 1.0 / FastMath.sqrt(3.0) + }; + + /** Weights for the 2 points method. */ + private static final double[] WEIGHTS_2 = { + 1.0, + 1.0 + }; + + /** Abscissas for the 3 points method. */ + private static final double[] ABSCISSAS_3 = { + -FastMath.sqrt(0.6), + 0.0, + FastMath.sqrt(0.6) + }; + + /** Weights for the 3 points method. */ + private static final double[] WEIGHTS_3 = { + 5.0 / 9.0, + 8.0 / 9.0, + 5.0 / 9.0 + }; + + /** Abscissas for the 4 points method. */ + private static final double[] ABSCISSAS_4 = { + -FastMath.sqrt((15.0 + 2.0 * FastMath.sqrt(30.0)) / 35.0), + -FastMath.sqrt((15.0 - 2.0 * FastMath.sqrt(30.0)) / 35.0), + FastMath.sqrt((15.0 - 2.0 * FastMath.sqrt(30.0)) / 35.0), + FastMath.sqrt((15.0 + 2.0 * FastMath.sqrt(30.0)) / 35.0) + }; + + /** Weights for the 4 points method. */ + private static final double[] WEIGHTS_4 = { + (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0, + (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0, + (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0, + (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0 + }; + + /** Abscissas for the 5 points method. */ + private static final double[] ABSCISSAS_5 = { + -FastMath.sqrt((35.0 + 2.0 * FastMath.sqrt(70.0)) / 63.0), + -FastMath.sqrt((35.0 - 2.0 * FastMath.sqrt(70.0)) / 63.0), + 0.0, + FastMath.sqrt((35.0 - 2.0 * FastMath.sqrt(70.0)) / 63.0), + FastMath.sqrt((35.0 + 2.0 * FastMath.sqrt(70.0)) / 63.0) + }; + + /** Weights for the 5 points method. */ + private static final double[] WEIGHTS_5 = { + (322.0 - 13.0 * FastMath.sqrt(70.0)) / 900.0, + (322.0 + 13.0 * FastMath.sqrt(70.0)) / 900.0, + 128.0 / 225.0, + (322.0 + 13.0 * FastMath.sqrt(70.0)) / 900.0, + (322.0 - 13.0 * FastMath.sqrt(70.0)) / 900.0 + }; + + /** Abscissas for the current method. */ + private final double[] abscissas; + + /** Weights for the current method. */ + private final double[] weights; + + /** + * Build a Legendre-Gauss integrator. + * @param n number of points desired (must be between 2 and 5 inclusive) + * @param defaultMaximalIterationCount maximum number of iterations + * @exception IllegalArgumentException if the number of points is not + * in the supported range + */ + public LegendreGaussIntegrator(final int n, final int defaultMaximalIterationCount) + throws IllegalArgumentException { + super(defaultMaximalIterationCount); + switch(n) { + case 2 : + abscissas = ABSCISSAS_2; + weights = WEIGHTS_2; + break; + case 3 : + abscissas = ABSCISSAS_3; + weights = WEIGHTS_3; + break; + case 4 : + abscissas = ABSCISSAS_4; + weights = WEIGHTS_4; + break; + case 5 : + abscissas = ABSCISSAS_5; + weights = WEIGHTS_5; + break; + default : + throw MathRuntimeException.createIllegalArgumentException( + LocalizedFormats.N_POINTS_GAUSS_LEGENDRE_INTEGRATOR_NOT_SUPPORTED, + n, 2, 5); + } + + } + + /** {@inheritDoc} */ + @Deprecated + public double integrate(final double min, final double max) + throws ConvergenceException, FunctionEvaluationException, IllegalArgumentException { + return integrate(f, min, max); + } + + /** {@inheritDoc} */ + public double integrate(final UnivariateRealFunction f, final double min, final double max) + throws ConvergenceException, FunctionEvaluationException, IllegalArgumentException { + + clearResult(); + verifyInterval(min, max); + verifyIterationCount(); + + // compute first estimate with a single step + double oldt = stage(f, min, max, 1); + + int n = 2; + for (int i = 0; i < maximalIterationCount; ++i) { + + // improve integral with a larger number of steps + final double t = stage(f, min, max, n); + + // estimate error + final double delta = FastMath.abs(t - oldt); + final double limit = + FastMath.max(absoluteAccuracy, + relativeAccuracy * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5); + + // check convergence + if ((i + 1 >= minimalIterationCount) && (delta <= limit)) { + setResult(t, i); + return result; + } + + // prepare next iteration + double ratio = FastMath.min(4, FastMath.pow(delta / limit, 0.5 / abscissas.length)); + n = FastMath.max((int) (ratio * n), n + 1); + oldt = t; + + } + + throw new MaxIterationsExceededException(maximalIterationCount); + + } + + /** + * Compute the n-th stage integral. + * @param f the integrand function + * @param min the lower bound for the interval + * @param max the upper bound for the interval + * @param n number of steps + * @return the value of n-th stage integral + * @throws FunctionEvaluationException if an error occurs evaluating the + * function + */ + private double stage(final UnivariateRealFunction f, + final double min, final double max, final int n) + throws FunctionEvaluationException { + + // set up the step for the current stage + final double step = (max - min) / n; + final double halfStep = step / 2.0; + + // integrate over all elementary steps + double midPoint = min + halfStep; + double sum = 0.0; + for (int i = 0; i < n; ++i) { + for (int j = 0; j < abscissas.length; ++j) { + sum += weights[j] * f.value(midPoint + halfStep * abscissas[j]); + } + midPoint += step; + } + + return halfStep * sum; + + } + +} |