Java tutorial
/* * 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 fvr. 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; } }