diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/analysis/integration')
17 files changed, 2576 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/BaseAbstractUnivariateIntegrator.java b/src/main/java/org/apache/commons/math3/analysis/integration/BaseAbstractUnivariateIntegrator.java new file mode 100644 index 0000000..74b959b --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/integration/BaseAbstractUnivariateIntegrator.java @@ -0,0 +1,297 @@ +/* + * 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.math3.analysis.integration; + +import org.apache.commons.math3.analysis.UnivariateFunction; +import org.apache.commons.math3.analysis.solvers.UnivariateSolverUtils; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.exception.TooManyEvaluationsException; +import org.apache.commons.math3.util.IntegerSequence; +import org.apache.commons.math3.util.MathUtils; + +/** + * Provide a default implementation for several generic functions. + * + * @since 1.2 + */ +public abstract class BaseAbstractUnivariateIntegrator implements UnivariateIntegrator { + + /** Default absolute accuracy. */ + public static final double DEFAULT_ABSOLUTE_ACCURACY = 1.0e-15; + + /** Default relative accuracy. */ + public static final double DEFAULT_RELATIVE_ACCURACY = 1.0e-6; + + /** Default minimal iteration count. */ + public static final int DEFAULT_MIN_ITERATIONS_COUNT = 3; + + /** Default maximal iteration count. */ + public static final int DEFAULT_MAX_ITERATIONS_COUNT = Integer.MAX_VALUE; + + /** The iteration count. + * @deprecated as of 3.6, this field has been replaced with {@link #incrementCount()} + */ + @Deprecated + protected org.apache.commons.math3.util.Incrementor iterations; + + /** The iteration count. */ + private IntegerSequence.Incrementor count; + + /** Maximum absolute error. */ + private final double absoluteAccuracy; + + /** Maximum relative error. */ + private final double relativeAccuracy; + + /** minimum number of iterations */ + private final int minimalIterationCount; + + /** The functions evaluation count. */ + private IntegerSequence.Incrementor evaluations; + + /** Function to integrate. */ + private UnivariateFunction function; + + /** Lower bound for the interval. */ + private double min; + + /** Upper bound for the interval. */ + private double max; + + /** + * Construct an integrator with given accuracies and iteration counts. + * <p> + * The meanings of the various parameters are: + * <ul> + * <li>relative accuracy: + * this is used to stop iterations if the absolute accuracy can't be + * achieved due to large values or short mantissa length. If this + * should be the primary criterion for convergence rather then a + * safety measure, set the absolute accuracy to a ridiculously small value, + * like {@link org.apache.commons.math3.util.Precision#SAFE_MIN Precision.SAFE_MIN}.</li> + * <li>absolute accuracy: + * The default is usually chosen so that results in the interval + * -10..-0.1 and +0.1..+10 can be found with a reasonable accuracy. If the + * expected absolute value of your results is of much smaller magnitude, set + * this to a smaller value.</li> + * <li>minimum number of iterations: + * minimal iteration is needed to avoid false early convergence, e.g. + * the sample points happen to be zeroes of the function. Users can + * use the default value or choose one that they see as appropriate.</li> + * <li>maximum number of iterations: + * usually a high iteration count indicates convergence problems. However, + * the "reasonable value" varies widely for different algorithms. Users are + * advised to use the default value supplied by the algorithm.</li> + * </ul> + * + * @param relativeAccuracy relative accuracy of the result + * @param absoluteAccuracy absolute accuracy of the result + * @param minimalIterationCount minimum number of iterations + * @param maximalIterationCount maximum number of iterations + * @exception NotStrictlyPositiveException if minimal number of iterations + * is not strictly positive + * @exception NumberIsTooSmallException if maximal number of iterations + * is lesser than or equal to the minimal number of iterations + */ + protected BaseAbstractUnivariateIntegrator(final double relativeAccuracy, + final double absoluteAccuracy, + final int minimalIterationCount, + final int maximalIterationCount) + throws NotStrictlyPositiveException, NumberIsTooSmallException { + + // accuracy settings + this.relativeAccuracy = relativeAccuracy; + this.absoluteAccuracy = absoluteAccuracy; + + // iterations count settings + if (minimalIterationCount <= 0) { + throw new NotStrictlyPositiveException(minimalIterationCount); + } + if (maximalIterationCount <= minimalIterationCount) { + throw new NumberIsTooSmallException(maximalIterationCount, minimalIterationCount, false); + } + this.minimalIterationCount = minimalIterationCount; + this.count = IntegerSequence.Incrementor.create().withMaximalCount(maximalIterationCount); + + @SuppressWarnings("deprecation") + org.apache.commons.math3.util.Incrementor wrapped = + org.apache.commons.math3.util.Incrementor.wrap(count); + this.iterations = wrapped; + + // prepare evaluations counter, but do not set it yet + evaluations = IntegerSequence.Incrementor.create(); + + } + + /** + * Construct an integrator with given accuracies. + * @param relativeAccuracy relative accuracy of the result + * @param absoluteAccuracy absolute accuracy of the result + */ + protected BaseAbstractUnivariateIntegrator(final double relativeAccuracy, + final double absoluteAccuracy) { + this(relativeAccuracy, absoluteAccuracy, + DEFAULT_MIN_ITERATIONS_COUNT, DEFAULT_MAX_ITERATIONS_COUNT); + } + + /** + * Construct an integrator with given iteration counts. + * @param minimalIterationCount minimum number of iterations + * @param maximalIterationCount maximum number of iterations + * @exception NotStrictlyPositiveException if minimal number of iterations + * is not strictly positive + * @exception NumberIsTooSmallException if maximal number of iterations + * is lesser than or equal to the minimal number of iterations + */ + protected BaseAbstractUnivariateIntegrator(final int minimalIterationCount, + final int maximalIterationCount) + throws NotStrictlyPositiveException, NumberIsTooSmallException { + this(DEFAULT_RELATIVE_ACCURACY, DEFAULT_ABSOLUTE_ACCURACY, + minimalIterationCount, maximalIterationCount); + } + + /** {@inheritDoc} */ + public double getRelativeAccuracy() { + return relativeAccuracy; + } + + /** {@inheritDoc} */ + public double getAbsoluteAccuracy() { + return absoluteAccuracy; + } + + /** {@inheritDoc} */ + public int getMinimalIterationCount() { + return minimalIterationCount; + } + + /** {@inheritDoc} */ + public int getMaximalIterationCount() { + return count.getMaximalCount(); + } + + /** {@inheritDoc} */ + public int getEvaluations() { + return evaluations.getCount(); + } + + /** {@inheritDoc} */ + public int getIterations() { + return count.getCount(); + } + + /** Increment the number of iterations. + * @exception MaxCountExceededException if the number of iterations + * exceeds the allowed maximum number + */ + protected void incrementCount() throws MaxCountExceededException { + count.increment(); + } + + /** + * @return the lower bound. + */ + protected double getMin() { + return min; + } + /** + * @return the upper bound. + */ + protected double getMax() { + return max; + } + + /** + * Compute the objective function value. + * + * @param point Point at which the objective function must be evaluated. + * @return the objective function value at specified point. + * @throws TooManyEvaluationsException if the maximal number of function + * evaluations is exceeded. + */ + protected double computeObjectiveValue(final double point) + throws TooManyEvaluationsException { + try { + evaluations.increment(); + } catch (MaxCountExceededException e) { + throw new TooManyEvaluationsException(e.getMax()); + } + return function.value(point); + } + + /** + * Prepare for computation. + * Subclasses must call this method if they override any of the + * {@code solve} methods. + * + * @param maxEval Maximum number of evaluations. + * @param f the integrand function + * @param lower the min bound for the interval + * @param upper the upper bound for the interval + * @throws NullArgumentException if {@code f} is {@code null}. + * @throws MathIllegalArgumentException if {@code min >= max}. + */ + protected void setup(final int maxEval, + final UnivariateFunction f, + final double lower, final double upper) + throws NullArgumentException, MathIllegalArgumentException { + + // Checks. + MathUtils.checkNotNull(f); + UnivariateSolverUtils.verifyInterval(lower, upper); + + // Reset. + min = lower; + max = upper; + function = f; + evaluations = evaluations.withMaximalCount(maxEval).withStart(0); + count = count.withStart(0); + + } + + /** {@inheritDoc} */ + public double integrate(final int maxEval, final UnivariateFunction f, + final double lower, final double upper) + throws TooManyEvaluationsException, MaxCountExceededException, + MathIllegalArgumentException, NullArgumentException { + + // Initialization. + setup(maxEval, f, lower, upper); + + // Perform computation. + return doIntegrate(); + + } + + /** + * Method for implementing actual integration algorithms in derived + * classes. + * + * @return the root. + * @throws TooManyEvaluationsException if the maximal number of evaluations + * is exceeded. + * @throws MaxCountExceededException if the maximum iteration count is exceeded + * or the integrator detects convergence problems otherwise + */ + protected abstract double doIntegrate() + throws TooManyEvaluationsException, MaxCountExceededException; + +} diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/IterativeLegendreGaussIntegrator.java b/src/main/java/org/apache/commons/math3/analysis/integration/IterativeLegendreGaussIntegrator.java new file mode 100644 index 0000000..20700dd --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/integration/IterativeLegendreGaussIntegrator.java @@ -0,0 +1,183 @@ +/* + * 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.math3.analysis.integration; + +import org.apache.commons.math3.analysis.UnivariateFunction; +import org.apache.commons.math3.analysis.integration.gauss.GaussIntegratorFactory; +import org.apache.commons.math3.analysis.integration.gauss.GaussIntegrator; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.exception.TooManyEvaluationsException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.util.FastMath; + +/** + * This algorithm divides the integration interval into equally-sized + * sub-interval and on each of them performs a + * <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html"> + * Legendre-Gauss</a> quadrature. + * Because of its <em>non-adaptive</em> nature, this algorithm can + * converge to a wrong value for the integral (for example, if the + * function is significantly different from zero toward the ends of the + * integration interval). + * In particular, a change of variables aimed at estimating integrals + * over infinite intervals as proposed + * <a href="http://en.wikipedia.org/w/index.php?title=Numerical_integration#Integrals_over_infinite_intervals"> + * here</a> should be avoided when using this class. + * + * @since 3.1 + */ + +public class IterativeLegendreGaussIntegrator + extends BaseAbstractUnivariateIntegrator { + /** Factory that computes the points and weights. */ + private static final GaussIntegratorFactory FACTORY + = new GaussIntegratorFactory(); + /** Number of integration points (per interval). */ + private final int numberOfPoints; + + /** + * Builds an integrator with given accuracies and iterations counts. + * + * @param n Number of integration points. + * @param relativeAccuracy Relative accuracy of the result. + * @param absoluteAccuracy Absolute accuracy of the result. + * @param minimalIterationCount Minimum number of iterations. + * @param maximalIterationCount Maximum number of iterations. + * @throws NotStrictlyPositiveException if minimal number of iterations + * or number of points are not strictly positive. + * @throws NumberIsTooSmallException if maximal number of iterations + * is smaller than or equal to the minimal number of iterations. + */ + public IterativeLegendreGaussIntegrator(final int n, + final double relativeAccuracy, + final double absoluteAccuracy, + final int minimalIterationCount, + final int maximalIterationCount) + throws NotStrictlyPositiveException, NumberIsTooSmallException { + super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount); + if (n <= 0) { + throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_POINTS, n); + } + numberOfPoints = n; + } + + /** + * Builds an integrator with given accuracies. + * + * @param n Number of integration points. + * @param relativeAccuracy Relative accuracy of the result. + * @param absoluteAccuracy Absolute accuracy of the result. + * @throws NotStrictlyPositiveException if {@code n < 1}. + */ + public IterativeLegendreGaussIntegrator(final int n, + final double relativeAccuracy, + final double absoluteAccuracy) + throws NotStrictlyPositiveException { + this(n, relativeAccuracy, absoluteAccuracy, + DEFAULT_MIN_ITERATIONS_COUNT, DEFAULT_MAX_ITERATIONS_COUNT); + } + + /** + * Builds an integrator with given iteration counts. + * + * @param n Number of integration points. + * @param minimalIterationCount Minimum number of iterations. + * @param maximalIterationCount Maximum number of iterations. + * @throws NotStrictlyPositiveException if minimal number of iterations + * is not strictly positive. + * @throws NumberIsTooSmallException if maximal number of iterations + * is smaller than or equal to the minimal number of iterations. + * @throws NotStrictlyPositiveException if {@code n < 1}. + */ + public IterativeLegendreGaussIntegrator(final int n, + final int minimalIterationCount, + final int maximalIterationCount) + throws NotStrictlyPositiveException, NumberIsTooSmallException { + this(n, DEFAULT_RELATIVE_ACCURACY, DEFAULT_ABSOLUTE_ACCURACY, + minimalIterationCount, maximalIterationCount); + } + + /** {@inheritDoc} */ + @Override + protected double doIntegrate() + throws MathIllegalArgumentException, TooManyEvaluationsException, MaxCountExceededException { + // Compute first estimate with a single step. + double oldt = stage(1); + + int n = 2; + while (true) { + // Improve integral with a larger number of steps. + final double t = stage(n); + + // Estimate the error. + final double delta = FastMath.abs(t - oldt); + final double limit = + FastMath.max(getAbsoluteAccuracy(), + getRelativeAccuracy() * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5); + + // check convergence + if (getIterations() + 1 >= getMinimalIterationCount() && + delta <= limit) { + return t; + } + + // Prepare next iteration. + final double ratio = FastMath.min(4, FastMath.pow(delta / limit, 0.5 / numberOfPoints)); + n = FastMath.max((int) (ratio * n), n + 1); + oldt = t; + incrementCount(); + } + } + + /** + * Compute the n-th stage integral. + * + * @param n Number of steps. + * @return the value of n-th stage integral. + * @throws TooManyEvaluationsException if the maximum number of evaluations + * is exceeded. + */ + private double stage(final int n) + throws TooManyEvaluationsException { + // Function to be integrated is stored in the base class. + final UnivariateFunction f = new UnivariateFunction() { + /** {@inheritDoc} */ + public double value(double x) + throws MathIllegalArgumentException, TooManyEvaluationsException { + return computeObjectiveValue(x); + } + }; + + final double min = getMin(); + final double max = getMax(); + final double step = (max - min) / n; + + double sum = 0; + for (int i = 0; i < n; i++) { + // Integrate over each sub-interval [a, b]. + final double a = min + i * step; + final double b = a + step; + final GaussIntegrator g = FACTORY.legendreHighPrecision(numberOfPoints, a, b); + sum += g.integrate(f); + } + + return sum; + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/LegendreGaussIntegrator.java b/src/main/java/org/apache/commons/math3/analysis/integration/LegendreGaussIntegrator.java new file mode 100644 index 0000000..bfb24a1 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/integration/LegendreGaussIntegrator.java @@ -0,0 +1,265 @@ +/* + * 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.math3.analysis.integration; + +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.exception.TooManyEvaluationsException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.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 function evaluations. A + * Legendre-Gauss integrator using an n-points quadrature formula can + * integrate 2n-1 degree polynomials exactly. + * </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> + * @since 1.2 + * @deprecated As of 3.1 (to be removed in 4.0). Please use + * {@link IterativeLegendreGaussIntegrator} instead. + */ +@Deprecated +public class LegendreGaussIntegrator extends BaseAbstractUnivariateIntegrator { + + /** 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 with given accuracies and iterations counts. + * @param n number of points desired (must be between 2 and 5 inclusive) + * @param relativeAccuracy relative accuracy of the result + * @param absoluteAccuracy absolute accuracy of the result + * @param minimalIterationCount minimum number of iterations + * @param maximalIterationCount maximum number of iterations + * @exception MathIllegalArgumentException if number of points is out of [2; 5] + * @exception NotStrictlyPositiveException if minimal number of iterations + * is not strictly positive + * @exception NumberIsTooSmallException if maximal number of iterations + * is lesser than or equal to the minimal number of iterations + */ + public LegendreGaussIntegrator(final int n, + final double relativeAccuracy, + final double absoluteAccuracy, + final int minimalIterationCount, + final int maximalIterationCount) + throws MathIllegalArgumentException, NotStrictlyPositiveException, NumberIsTooSmallException { + super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount); + 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 new MathIllegalArgumentException( + LocalizedFormats.N_POINTS_GAUSS_LEGENDRE_INTEGRATOR_NOT_SUPPORTED, + n, 2, 5); + } + + } + + /** + * Build a Legendre-Gauss integrator with given accuracies. + * @param n number of points desired (must be between 2 and 5 inclusive) + * @param relativeAccuracy relative accuracy of the result + * @param absoluteAccuracy absolute accuracy of the result + * @exception MathIllegalArgumentException if number of points is out of [2; 5] + */ + public LegendreGaussIntegrator(final int n, + final double relativeAccuracy, + final double absoluteAccuracy) + throws MathIllegalArgumentException { + this(n, relativeAccuracy, absoluteAccuracy, + DEFAULT_MIN_ITERATIONS_COUNT, DEFAULT_MAX_ITERATIONS_COUNT); + } + + /** + * Build a Legendre-Gauss integrator with given iteration counts. + * @param n number of points desired (must be between 2 and 5 inclusive) + * @param minimalIterationCount minimum number of iterations + * @param maximalIterationCount maximum number of iterations + * @exception MathIllegalArgumentException if number of points is out of [2; 5] + * @exception NotStrictlyPositiveException if minimal number of iterations + * is not strictly positive + * @exception NumberIsTooSmallException if maximal number of iterations + * is lesser than or equal to the minimal number of iterations + */ + public LegendreGaussIntegrator(final int n, + final int minimalIterationCount, + final int maximalIterationCount) + throws MathIllegalArgumentException { + this(n, DEFAULT_RELATIVE_ACCURACY, DEFAULT_ABSOLUTE_ACCURACY, + minimalIterationCount, maximalIterationCount); + } + + /** {@inheritDoc} */ + @Override + protected double doIntegrate() + throws MathIllegalArgumentException, TooManyEvaluationsException, MaxCountExceededException { + + // compute first estimate with a single step + double oldt = stage(1); + + int n = 2; + while (true) { + + // improve integral with a larger number of steps + final double t = stage(n); + + // estimate error + final double delta = FastMath.abs(t - oldt); + final double limit = + FastMath.max(getAbsoluteAccuracy(), + getRelativeAccuracy() * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5); + + // check convergence + if ((getIterations() + 1 >= getMinimalIterationCount()) && (delta <= limit)) { + return t; + } + + // 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; + incrementCount(); + + } + + } + + /** + * Compute the n-th stage integral. + * @param n number of steps + * @return the value of n-th stage integral + * @throws TooManyEvaluationsException if the maximum number of evaluations + * is exceeded. + */ + private double stage(final int n) + throws TooManyEvaluationsException { + + // set up the step for the current stage + final double step = (getMax() - getMin()) / n; + final double halfStep = step / 2.0; + + // integrate over all elementary steps + double midPoint = getMin() + halfStep; + double sum = 0.0; + for (int i = 0; i < n; ++i) { + for (int j = 0; j < abscissas.length; ++j) { + sum += weights[j] * computeObjectiveValue(midPoint + halfStep * abscissas[j]); + } + midPoint += step; + } + + return halfStep * sum; + + } + +} diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/MidPointIntegrator.java b/src/main/java/org/apache/commons/math3/analysis/integration/MidPointIntegrator.java new file mode 100644 index 0000000..766a917 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/integration/MidPointIntegrator.java @@ -0,0 +1,169 @@ +/* + * 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.math3.analysis.integration; + +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.NumberIsTooLargeException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.exception.TooManyEvaluationsException; +import org.apache.commons.math3.util.FastMath; + +/** + * Implements the <a href="http://en.wikipedia.org/wiki/Midpoint_method"> + * Midpoint Rule</a> for integration of real univariate functions. For + * reference, see <b>Numerical Mathematics</b>, ISBN 0387989595, + * chapter 9.2. + * <p> + * The function should be integrable.</p> + * + * @since 3.3 + */ +public class MidPointIntegrator extends BaseAbstractUnivariateIntegrator { + + /** Maximum number of iterations for midpoint. */ + public static final int MIDPOINT_MAX_ITERATIONS_COUNT = 64; + + /** + * Build a midpoint integrator with given accuracies and iterations counts. + * @param relativeAccuracy relative accuracy of the result + * @param absoluteAccuracy absolute accuracy of the result + * @param minimalIterationCount minimum number of iterations + * @param maximalIterationCount maximum number of iterations + * (must be less than or equal to {@link #MIDPOINT_MAX_ITERATIONS_COUNT} + * @exception NotStrictlyPositiveException if minimal number of iterations + * is not strictly positive + * @exception NumberIsTooSmallException if maximal number of iterations + * is lesser than or equal to the minimal number of iterations + * @exception NumberIsTooLargeException if maximal number of iterations + * is greater than {@link #MIDPOINT_MAX_ITERATIONS_COUNT} + */ + public MidPointIntegrator(final double relativeAccuracy, + final double absoluteAccuracy, + final int minimalIterationCount, + final int maximalIterationCount) + throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException { + super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount); + if (maximalIterationCount > MIDPOINT_MAX_ITERATIONS_COUNT) { + throw new NumberIsTooLargeException(maximalIterationCount, + MIDPOINT_MAX_ITERATIONS_COUNT, false); + } + } + + /** + * Build a midpoint integrator with given iteration counts. + * @param minimalIterationCount minimum number of iterations + * @param maximalIterationCount maximum number of iterations + * (must be less than or equal to {@link #MIDPOINT_MAX_ITERATIONS_COUNT} + * @exception NotStrictlyPositiveException if minimal number of iterations + * is not strictly positive + * @exception NumberIsTooSmallException if maximal number of iterations + * is lesser than or equal to the minimal number of iterations + * @exception NumberIsTooLargeException if maximal number of iterations + * is greater than {@link #MIDPOINT_MAX_ITERATIONS_COUNT} + */ + public MidPointIntegrator(final int minimalIterationCount, + final int maximalIterationCount) + throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException { + super(minimalIterationCount, maximalIterationCount); + if (maximalIterationCount > MIDPOINT_MAX_ITERATIONS_COUNT) { + throw new NumberIsTooLargeException(maximalIterationCount, + MIDPOINT_MAX_ITERATIONS_COUNT, false); + } + } + + /** + * Construct a midpoint integrator with default settings. + * (max iteration count set to {@link #MIDPOINT_MAX_ITERATIONS_COUNT}) + */ + public MidPointIntegrator() { + super(DEFAULT_MIN_ITERATIONS_COUNT, MIDPOINT_MAX_ITERATIONS_COUNT); + } + + /** + * Compute the n-th stage integral of midpoint rule. + * This function should only be called by API <code>integrate()</code> in the package. + * To save time it does not verify arguments - caller does. + * <p> + * The interval is divided equally into 2^n sections rather than an + * arbitrary m sections because this configuration can best utilize the + * already computed values.</p> + * + * @param n the stage of 1/2 refinement. Must be larger than 0. + * @param previousStageResult Result from the previous call to the + * {@code stage} method. + * @param min Lower bound of the integration interval. + * @param diffMaxMin Difference between the lower bound and upper bound + * of the integration interval. + * @return the value of n-th stage integral + * @throws TooManyEvaluationsException if the maximal number of evaluations + * is exceeded. + */ + private double stage(final int n, + double previousStageResult, + double min, + double diffMaxMin) + throws TooManyEvaluationsException { + + // number of new points in this stage + final long np = 1L << (n - 1); + double sum = 0; + + // spacing between adjacent new points + final double spacing = diffMaxMin / np; + + // the first new point + double x = min + 0.5 * spacing; + for (long i = 0; i < np; i++) { + sum += computeObjectiveValue(x); + x += spacing; + } + // add the new sum to previously calculated result + return 0.5 * (previousStageResult + sum * spacing); + } + + + /** {@inheritDoc} */ + @Override + protected double doIntegrate() + throws MathIllegalArgumentException, TooManyEvaluationsException, MaxCountExceededException { + + final double min = getMin(); + final double diff = getMax() - min; + final double midPoint = min + 0.5 * diff; + + double oldt = diff * computeObjectiveValue(midPoint); + + while (true) { + incrementCount(); + final int i = getIterations(); + final double t = stage(i, oldt, min, diff); + if (i >= getMinimalIterationCount()) { + final double delta = FastMath.abs(t - oldt); + final double rLimit = + getRelativeAccuracy() * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5; + if ((delta <= rLimit) || (delta <= getAbsoluteAccuracy())) { + return t; + } + } + oldt = t; + } + + } + +} diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/RombergIntegrator.java b/src/main/java/org/apache/commons/math3/analysis/integration/RombergIntegrator.java new file mode 100644 index 0000000..125d251 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/integration/RombergIntegrator.java @@ -0,0 +1,142 @@ +/* + * 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.math3.analysis.integration; + +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.NumberIsTooLargeException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.exception.TooManyEvaluationsException; +import org.apache.commons.math3.util.FastMath; + +/** + * Implements the <a href="http://mathworld.wolfram.com/RombergIntegration.html"> + * Romberg Algorithm</a> for integration of real univariate functions. For + * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X, + * chapter 3. + * <p> + * Romberg integration employs k successive refinements of the trapezoid + * rule to remove error terms less than order O(N^(-2k)). Simpson's rule + * is a special case of k = 2.</p> + * + * @since 1.2 + */ +public class RombergIntegrator extends BaseAbstractUnivariateIntegrator { + + /** Maximal number of iterations for Romberg. */ + public static final int ROMBERG_MAX_ITERATIONS_COUNT = 32; + + /** + * Build a Romberg integrator with given accuracies and iterations counts. + * @param relativeAccuracy relative accuracy of the result + * @param absoluteAccuracy absolute accuracy of the result + * @param minimalIterationCount minimum number of iterations + * @param maximalIterationCount maximum number of iterations + * (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT}) + * @exception NotStrictlyPositiveException if minimal number of iterations + * is not strictly positive + * @exception NumberIsTooSmallException if maximal number of iterations + * is lesser than or equal to the minimal number of iterations + * @exception NumberIsTooLargeException if maximal number of iterations + * is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT} + */ + public RombergIntegrator(final double relativeAccuracy, + final double absoluteAccuracy, + final int minimalIterationCount, + final int maximalIterationCount) + throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException { + super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount); + if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) { + throw new NumberIsTooLargeException(maximalIterationCount, + ROMBERG_MAX_ITERATIONS_COUNT, false); + } + } + + /** + * Build a Romberg integrator with given iteration counts. + * @param minimalIterationCount minimum number of iterations + * @param maximalIterationCount maximum number of iterations + * (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT}) + * @exception NotStrictlyPositiveException if minimal number of iterations + * is not strictly positive + * @exception NumberIsTooSmallException if maximal number of iterations + * is lesser than or equal to the minimal number of iterations + * @exception NumberIsTooLargeException if maximal number of iterations + * is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT} + */ + public RombergIntegrator(final int minimalIterationCount, + final int maximalIterationCount) + throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException { + super(minimalIterationCount, maximalIterationCount); + if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) { + throw new NumberIsTooLargeException(maximalIterationCount, + ROMBERG_MAX_ITERATIONS_COUNT, false); + } + } + + /** + * Construct a Romberg integrator with default settings + * (max iteration count set to {@link #ROMBERG_MAX_ITERATIONS_COUNT}) + */ + public RombergIntegrator() { + super(DEFAULT_MIN_ITERATIONS_COUNT, ROMBERG_MAX_ITERATIONS_COUNT); + } + + /** {@inheritDoc} */ + @Override + protected double doIntegrate() + throws TooManyEvaluationsException, MaxCountExceededException { + + final int m = getMaximalIterationCount() + 1; + double previousRow[] = new double[m]; + double currentRow[] = new double[m]; + + TrapezoidIntegrator qtrap = new TrapezoidIntegrator(); + currentRow[0] = qtrap.stage(this, 0); + incrementCount(); + double olds = currentRow[0]; + while (true) { + + final int i = getIterations(); + + // switch rows + final double[] tmpRow = previousRow; + previousRow = currentRow; + currentRow = tmpRow; + + currentRow[0] = qtrap.stage(this, i); + incrementCount(); + for (int j = 1; j <= i; j++) { + // Richardson extrapolation coefficient + final double r = (1L << (2 * j)) - 1; + final double tIJm1 = currentRow[j - 1]; + currentRow[j] = tIJm1 + (tIJm1 - previousRow[j - 1]) / r; + } + final double s = currentRow[i]; + if (i >= getMinimalIterationCount()) { + final double delta = FastMath.abs(s - olds); + final double rLimit = getRelativeAccuracy() * (FastMath.abs(olds) + FastMath.abs(s)) * 0.5; + if ((delta <= rLimit) || (delta <= getAbsoluteAccuracy())) { + return s; + } + } + olds = s; + } + + } + +} diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/SimpsonIntegrator.java b/src/main/java/org/apache/commons/math3/analysis/integration/SimpsonIntegrator.java new file mode 100644 index 0000000..527bb82 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/integration/SimpsonIntegrator.java @@ -0,0 +1,129 @@ +/* + * 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.math3.analysis.integration; + +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.NumberIsTooLargeException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.exception.TooManyEvaluationsException; +import org.apache.commons.math3.util.FastMath; + +/** + * Implements <a href="http://mathworld.wolfram.com/SimpsonsRule.html"> + * Simpson's Rule</a> for integration of real univariate functions. For + * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X, + * chapter 3. + * <p> + * This implementation employs the basic trapezoid rule to calculate Simpson's + * rule.</p> + * + * @since 1.2 + */ +public class SimpsonIntegrator extends BaseAbstractUnivariateIntegrator { + + /** Maximal number of iterations for Simpson. */ + public static final int SIMPSON_MAX_ITERATIONS_COUNT = 64; + + /** + * Build a Simpson integrator with given accuracies and iterations counts. + * @param relativeAccuracy relative accuracy of the result + * @param absoluteAccuracy absolute accuracy of the result + * @param minimalIterationCount minimum number of iterations + * @param maximalIterationCount maximum number of iterations + * (must be less than or equal to {@link #SIMPSON_MAX_ITERATIONS_COUNT}) + * @exception NotStrictlyPositiveException if minimal number of iterations + * is not strictly positive + * @exception NumberIsTooSmallException if maximal number of iterations + * is lesser than or equal to the minimal number of iterations + * @exception NumberIsTooLargeException if maximal number of iterations + * is greater than {@link #SIMPSON_MAX_ITERATIONS_COUNT} + */ + public SimpsonIntegrator(final double relativeAccuracy, + final double absoluteAccuracy, + final int minimalIterationCount, + final int maximalIterationCount) + throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException { + super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount); + if (maximalIterationCount > SIMPSON_MAX_ITERATIONS_COUNT) { + throw new NumberIsTooLargeException(maximalIterationCount, + SIMPSON_MAX_ITERATIONS_COUNT, false); + } + } + + /** + * Build a Simpson integrator with given iteration counts. + * @param minimalIterationCount minimum number of iterations + * @param maximalIterationCount maximum number of iterations + * (must be less than or equal to {@link #SIMPSON_MAX_ITERATIONS_COUNT}) + * @exception NotStrictlyPositiveException if minimal number of iterations + * is not strictly positive + * @exception NumberIsTooSmallException if maximal number of iterations + * is lesser than or equal to the minimal number of iterations + * @exception NumberIsTooLargeException if maximal number of iterations + * is greater than {@link #SIMPSON_MAX_ITERATIONS_COUNT} + */ + public SimpsonIntegrator(final int minimalIterationCount, + final int maximalIterationCount) + throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException { + super(minimalIterationCount, maximalIterationCount); + if (maximalIterationCount > SIMPSON_MAX_ITERATIONS_COUNT) { + throw new NumberIsTooLargeException(maximalIterationCount, + SIMPSON_MAX_ITERATIONS_COUNT, false); + } + } + + /** + * Construct an integrator with default settings. + * (max iteration count set to {@link #SIMPSON_MAX_ITERATIONS_COUNT}) + */ + public SimpsonIntegrator() { + super(DEFAULT_MIN_ITERATIONS_COUNT, SIMPSON_MAX_ITERATIONS_COUNT); + } + + /** {@inheritDoc} */ + @Override + protected double doIntegrate() + throws TooManyEvaluationsException, MaxCountExceededException { + + TrapezoidIntegrator qtrap = new TrapezoidIntegrator(); + if (getMinimalIterationCount() == 1) { + return (4 * qtrap.stage(this, 1) - qtrap.stage(this, 0)) / 3.0; + } + + // Simpson's rule requires at least two trapezoid stages. + double olds = 0; + double oldt = qtrap.stage(this, 0); + while (true) { + final double t = qtrap.stage(this, getIterations()); + incrementCount(); + final double s = (4 * t - oldt) / 3.0; + if (getIterations() >= getMinimalIterationCount()) { + final double delta = FastMath.abs(s - olds); + final double rLimit = + getRelativeAccuracy() * (FastMath.abs(olds) + FastMath.abs(s)) * 0.5; + if ((delta <= rLimit) || (delta <= getAbsoluteAccuracy())) { + return s; + } + } + olds = s; + oldt = t; + } + + } + +} diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/TrapezoidIntegrator.java b/src/main/java/org/apache/commons/math3/analysis/integration/TrapezoidIntegrator.java new file mode 100644 index 0000000..8a737f1 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/integration/TrapezoidIntegrator.java @@ -0,0 +1,168 @@ +/* + * 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.math3.analysis.integration; + +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.NumberIsTooLargeException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.exception.TooManyEvaluationsException; +import org.apache.commons.math3.util.FastMath; + +/** + * Implements the <a href="http://mathworld.wolfram.com/TrapezoidalRule.html"> + * Trapezoid Rule</a> for integration of real univariate functions. For + * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X, + * chapter 3. + * <p> + * The function should be integrable.</p> + * + * @since 1.2 + */ +public class TrapezoidIntegrator extends BaseAbstractUnivariateIntegrator { + + /** Maximum number of iterations for trapezoid. */ + public static final int TRAPEZOID_MAX_ITERATIONS_COUNT = 64; + + /** Intermediate result. */ + private double s; + + /** + * Build a trapezoid integrator with given accuracies and iterations counts. + * @param relativeAccuracy relative accuracy of the result + * @param absoluteAccuracy absolute accuracy of the result + * @param minimalIterationCount minimum number of iterations + * @param maximalIterationCount maximum number of iterations + * (must be less than or equal to {@link #TRAPEZOID_MAX_ITERATIONS_COUNT} + * @exception NotStrictlyPositiveException if minimal number of iterations + * is not strictly positive + * @exception NumberIsTooSmallException if maximal number of iterations + * is lesser than or equal to the minimal number of iterations + * @exception NumberIsTooLargeException if maximal number of iterations + * is greater than {@link #TRAPEZOID_MAX_ITERATIONS_COUNT} + */ + public TrapezoidIntegrator(final double relativeAccuracy, + final double absoluteAccuracy, + final int minimalIterationCount, + final int maximalIterationCount) + throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException { + super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount); + if (maximalIterationCount > TRAPEZOID_MAX_ITERATIONS_COUNT) { + throw new NumberIsTooLargeException(maximalIterationCount, + TRAPEZOID_MAX_ITERATIONS_COUNT, false); + } + } + + /** + * Build a trapezoid integrator with given iteration counts. + * @param minimalIterationCount minimum number of iterations + * @param maximalIterationCount maximum number of iterations + * (must be less than or equal to {@link #TRAPEZOID_MAX_ITERATIONS_COUNT} + * @exception NotStrictlyPositiveException if minimal number of iterations + * is not strictly positive + * @exception NumberIsTooSmallException if maximal number of iterations + * is lesser than or equal to the minimal number of iterations + * @exception NumberIsTooLargeException if maximal number of iterations + * is greater than {@link #TRAPEZOID_MAX_ITERATIONS_COUNT} + */ + public TrapezoidIntegrator(final int minimalIterationCount, + final int maximalIterationCount) + throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException { + super(minimalIterationCount, maximalIterationCount); + if (maximalIterationCount > TRAPEZOID_MAX_ITERATIONS_COUNT) { + throw new NumberIsTooLargeException(maximalIterationCount, + TRAPEZOID_MAX_ITERATIONS_COUNT, false); + } + } + + /** + * Construct a trapezoid integrator with default settings. + * (max iteration count set to {@link #TRAPEZOID_MAX_ITERATIONS_COUNT}) + */ + public TrapezoidIntegrator() { + super(DEFAULT_MIN_ITERATIONS_COUNT, TRAPEZOID_MAX_ITERATIONS_COUNT); + } + + /** + * Compute the n-th stage integral of trapezoid rule. This function + * should only be called by API <code>integrate()</code> in the package. + * To save time it does not verify arguments - caller does. + * <p> + * The interval is divided equally into 2^n sections rather than an + * arbitrary m sections because this configuration can best utilize the + * already computed values.</p> + * + * @param baseIntegrator integrator holding integration parameters + * @param n the stage of 1/2 refinement, n = 0 is no refinement + * @return the value of n-th stage integral + * @throws TooManyEvaluationsException if the maximal number of evaluations + * is exceeded. + */ + double stage(final BaseAbstractUnivariateIntegrator baseIntegrator, final int n) + throws TooManyEvaluationsException { + + if (n == 0) { + final double max = baseIntegrator.getMax(); + final double min = baseIntegrator.getMin(); + s = 0.5 * (max - min) * + (baseIntegrator.computeObjectiveValue(min) + + baseIntegrator.computeObjectiveValue(max)); + return s; + } else { + final long np = 1L << (n-1); // number of new points in this stage + double sum = 0; + final double max = baseIntegrator.getMax(); + final double min = baseIntegrator.getMin(); + // spacing between adjacent new points + final double spacing = (max - min) / np; + double x = min + 0.5 * spacing; // the first new point + for (long i = 0; i < np; i++) { + sum += baseIntegrator.computeObjectiveValue(x); + x += spacing; + } + // add the new sum to previously calculated result + s = 0.5 * (s + sum * spacing); + return s; + } + } + + /** {@inheritDoc} */ + @Override + protected double doIntegrate() + throws MathIllegalArgumentException, TooManyEvaluationsException, MaxCountExceededException { + + double oldt = stage(this, 0); + incrementCount(); + while (true) { + final int i = getIterations(); + final double t = stage(this, i); + if (i >= getMinimalIterationCount()) { + final double delta = FastMath.abs(t - oldt); + final double rLimit = + getRelativeAccuracy() * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5; + if ((delta <= rLimit) || (delta <= getAbsoluteAccuracy())) { + return t; + } + } + oldt = t; + incrementCount(); + } + + } + +} diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/UnivariateIntegrator.java b/src/main/java/org/apache/commons/math3/analysis/integration/UnivariateIntegrator.java new file mode 100644 index 0000000..6453eba --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/integration/UnivariateIntegrator.java @@ -0,0 +1,95 @@ +/* + * 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.math3.analysis.integration; + +import org.apache.commons.math3.analysis.UnivariateFunction; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.exception.TooManyEvaluationsException; + +/** + * Interface for univariate real integration algorithms. + * + * @since 1.2 + */ +public interface UnivariateIntegrator { + + /** + * Get the relative accuracy. + * + * @return the accuracy + */ + double getRelativeAccuracy(); + + /** + * Get the absolute accuracy. + * + * @return the accuracy + */ + double getAbsoluteAccuracy(); + + /** + * Get the min limit for the number of iterations. + * + * @return the actual min limit + */ + int getMinimalIterationCount(); + + /** + * Get the upper limit for the number of iterations. + * + * @return the actual upper limit + */ + int getMaximalIterationCount(); + + /** + * Integrate the function in the given interval. + * + * @param maxEval Maximum number of evaluations. + * @param f the integrand function + * @param min the lower bound for the interval + * @param max the upper bound for the interval + * @return the value of integral + * @throws TooManyEvaluationsException if the maximum number of function + * evaluations is exceeded + * @throws MaxCountExceededException if the maximum iteration count is exceeded + * or the integrator detects convergence problems otherwise + * @throws MathIllegalArgumentException if {@code min > max} or the endpoints do not + * satisfy the requirements specified by the integrator + * @throws NullArgumentException if {@code f} is {@code null}. + */ + double integrate(int maxEval, UnivariateFunction f, double min, + double max) + throws TooManyEvaluationsException, MaxCountExceededException, + MathIllegalArgumentException, NullArgumentException; + + /** + * Get the number of function evaluations of the last run of the integrator. + * + * @return number of function evaluations + */ + int getEvaluations(); + + /** + * Get the number of iterations of the last run of the integrator. + * + * @return number of iterations + */ + int getIterations(); + +} diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/BaseRuleFactory.java b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/BaseRuleFactory.java new file mode 100644 index 0000000..556fa0c --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/BaseRuleFactory.java @@ -0,0 +1,153 @@ +/* + * 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.math3.analysis.integration.gauss; + +import java.util.Map; +import java.util.TreeMap; +import org.apache.commons.math3.util.Pair; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.util.LocalizedFormats; + +/** + * Base class for rules that determines the integration nodes and their + * weights. + * Subclasses must implement the {@link #computeRule(int) computeRule} method. + * + * @param <T> Type of the number used to represent the points and weights of + * the quadrature rules. + * + * @since 3.1 + */ +public abstract class BaseRuleFactory<T extends Number> { + /** List of points and weights, indexed by the order of the rule. */ + private final Map<Integer, Pair<T[], T[]>> pointsAndWeights + = new TreeMap<Integer, Pair<T[], T[]>>(); + /** Cache for double-precision rules. */ + private final Map<Integer, Pair<double[], double[]>> pointsAndWeightsDouble + = new TreeMap<Integer, Pair<double[], double[]>>(); + + /** + * Gets a copy of the quadrature rule with the given number of integration + * points. + * + * @param numberOfPoints Number of integration points. + * @return a copy of the integration rule. + * @throws NotStrictlyPositiveException if {@code numberOfPoints < 1}. + * @throws DimensionMismatchException if the elements of the rule pair do not + * have the same length. + */ + public Pair<double[], double[]> getRule(int numberOfPoints) + throws NotStrictlyPositiveException, DimensionMismatchException { + + if (numberOfPoints <= 0) { + throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_POINTS, + numberOfPoints); + } + + // Try to obtain the rule from the cache. + Pair<double[], double[]> cached = pointsAndWeightsDouble.get(numberOfPoints); + + if (cached == null) { + // Rule not computed yet. + + // Compute the rule. + final Pair<T[], T[]> rule = getRuleInternal(numberOfPoints); + cached = convertToDouble(rule); + + // Cache it. + pointsAndWeightsDouble.put(numberOfPoints, cached); + } + + // Return a copy. + return new Pair<double[], double[]>(cached.getFirst().clone(), + cached.getSecond().clone()); + } + + /** + * Gets a rule. + * Synchronization ensures that rules will be computed and added to the + * cache at most once. + * The returned rule is a reference into the cache. + * + * @param numberOfPoints Order of the rule to be retrieved. + * @return the points and weights corresponding to the given order. + * @throws DimensionMismatchException if the elements of the rule pair do not + * have the same length. + */ + protected synchronized Pair<T[], T[]> getRuleInternal(int numberOfPoints) + throws DimensionMismatchException { + final Pair<T[], T[]> rule = pointsAndWeights.get(numberOfPoints); + if (rule == null) { + addRule(computeRule(numberOfPoints)); + // The rule should be available now. + return getRuleInternal(numberOfPoints); + } + return rule; + } + + /** + * Stores a rule. + * + * @param rule Rule to be stored. + * @throws DimensionMismatchException if the elements of the pair do not + * have the same length. + */ + protected void addRule(Pair<T[], T[]> rule) throws DimensionMismatchException { + if (rule.getFirst().length != rule.getSecond().length) { + throw new DimensionMismatchException(rule.getFirst().length, + rule.getSecond().length); + } + + pointsAndWeights.put(rule.getFirst().length, rule); + } + + /** + * Computes the rule for the given order. + * + * @param numberOfPoints Order of the rule to be computed. + * @return the computed rule. + * @throws DimensionMismatchException if the elements of the pair do not + * have the same length. + */ + protected abstract Pair<T[], T[]> computeRule(int numberOfPoints) + throws DimensionMismatchException; + + /** + * Converts the from the actual {@code Number} type to {@code double} + * + * @param <T> Type of the number used to represent the points and + * weights of the quadrature rules. + * @param rule Points and weights. + * @return points and weights as {@code double}s. + */ + private static <T extends Number> Pair<double[], double[]> convertToDouble(Pair<T[], T[]> rule) { + final T[] pT = rule.getFirst(); + final T[] wT = rule.getSecond(); + + final int len = pT.length; + final double[] pD = new double[len]; + final double[] wD = new double[len]; + + for (int i = 0; i < len; i++) { + pD[i] = pT[i].doubleValue(); + wD[i] = wT[i].doubleValue(); + } + + return new Pair<double[], double[]>(pD, wD); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegrator.java b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegrator.java new file mode 100644 index 0000000..5c7b37f --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegrator.java @@ -0,0 +1,129 @@ +/* + * 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.math3.analysis.integration.gauss; + +import org.apache.commons.math3.analysis.UnivariateFunction; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NonMonotonicSequenceException; +import org.apache.commons.math3.util.MathArrays; +import org.apache.commons.math3.util.Pair; + +/** + * Class that implements the Gaussian rule for + * {@link #integrate(UnivariateFunction) integrating} a weighted + * function. + * + * @since 3.1 + */ +public class GaussIntegrator { + /** Nodes. */ + private final double[] points; + /** Nodes weights. */ + private final double[] weights; + + /** + * Creates an integrator from the given {@code points} and {@code weights}. + * The integration interval is defined by the first and last value of + * {@code points} which must be sorted in increasing order. + * + * @param points Integration points. + * @param weights Weights of the corresponding integration nodes. + * @throws NonMonotonicSequenceException if the {@code points} are not + * sorted in increasing order. + * @throws DimensionMismatchException if points and weights don't have the same length + */ + public GaussIntegrator(double[] points, + double[] weights) + throws NonMonotonicSequenceException, DimensionMismatchException { + if (points.length != weights.length) { + throw new DimensionMismatchException(points.length, + weights.length); + } + + MathArrays.checkOrder(points, MathArrays.OrderDirection.INCREASING, true, true); + + this.points = points.clone(); + this.weights = weights.clone(); + } + + /** + * Creates an integrator from the given pair of points (first element of + * the pair) and weights (second element of the pair. + * + * @param pointsAndWeights Integration points and corresponding weights. + * @throws NonMonotonicSequenceException if the {@code points} are not + * sorted in increasing order. + * + * @see #GaussIntegrator(double[], double[]) + */ + public GaussIntegrator(Pair<double[], double[]> pointsAndWeights) + throws NonMonotonicSequenceException { + this(pointsAndWeights.getFirst(), pointsAndWeights.getSecond()); + } + + /** + * Returns an estimate of the integral of {@code f(x) * w(x)}, + * where {@code w} is a weight function that depends on the actual + * flavor of the Gauss integration scheme. + * The algorithm uses the points and associated weights, as passed + * to the {@link #GaussIntegrator(double[],double[]) constructor}. + * + * @param f Function to integrate. + * @return the integral of the weighted function. + */ + public double integrate(UnivariateFunction f) { + double s = 0; + double c = 0; + for (int i = 0; i < points.length; i++) { + final double x = points[i]; + final double w = weights[i]; + final double y = w * f.value(x) - c; + final double t = s + y; + c = (t - s) - y; + s = t; + } + return s; + } + + /** + * @return the order of the integration rule (the number of integration + * points). + */ + public int getNumberOfPoints() { + return points.length; + } + + /** + * Gets the integration point at the given index. + * The index must be in the valid range but no check is performed. + * @param index index of the integration point + * @return the integration point. + */ + public double getPoint(int index) { + return points[index]; + } + + /** + * Gets the weight of the integration point at the given index. + * The index must be in the valid range but no check is performed. + * @param index index of the integration point + * @return the weight. + */ + public double getWeight(int index) { + return weights[index]; + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegratorFactory.java b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegratorFactory.java new file mode 100644 index 0000000..570a7ba --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegratorFactory.java @@ -0,0 +1,167 @@ +/* + * 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.math3.analysis.integration.gauss; + +import java.math.BigDecimal; + +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.util.Pair; + +/** + * Class that provides different ways to compute the nodes and weights to be + * used by the {@link GaussIntegrator Gaussian integration rule}. + * + * @since 3.1 + */ +public class GaussIntegratorFactory { + /** Generator of Gauss-Legendre integrators. */ + private final BaseRuleFactory<Double> legendre = new LegendreRuleFactory(); + /** Generator of Gauss-Legendre integrators. */ + private final BaseRuleFactory<BigDecimal> legendreHighPrecision = new LegendreHighPrecisionRuleFactory(); + /** Generator of Gauss-Hermite integrators. */ + private final BaseRuleFactory<Double> hermite = new HermiteRuleFactory(); + + /** + * Creates a Gauss-Legendre integrator of the given order. + * The call to the + * {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction) + * integrate} method will perform an integration on the natural interval + * {@code [-1 , 1]}. + * + * @param numberOfPoints Order of the integration rule. + * @return a Gauss-Legendre integrator. + */ + public GaussIntegrator legendre(int numberOfPoints) { + return new GaussIntegrator(getRule(legendre, numberOfPoints)); + } + + /** + * Creates a Gauss-Legendre integrator of the given order. + * The call to the + * {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction) + * integrate} method will perform an integration on the given interval. + * + * @param numberOfPoints Order of the integration rule. + * @param lowerBound Lower bound of the integration interval. + * @param upperBound Upper bound of the integration interval. + * @return a Gauss-Legendre integrator. + * @throws NotStrictlyPositiveException if number of points is not positive + */ + public GaussIntegrator legendre(int numberOfPoints, + double lowerBound, + double upperBound) + throws NotStrictlyPositiveException { + return new GaussIntegrator(transform(getRule(legendre, numberOfPoints), + lowerBound, upperBound)); + } + + /** + * Creates a Gauss-Legendre integrator of the given order. + * The call to the + * {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction) + * integrate} method will perform an integration on the natural interval + * {@code [-1 , 1]}. + * + * @param numberOfPoints Order of the integration rule. + * @return a Gauss-Legendre integrator. + * @throws NotStrictlyPositiveException if number of points is not positive + */ + public GaussIntegrator legendreHighPrecision(int numberOfPoints) + throws NotStrictlyPositiveException { + return new GaussIntegrator(getRule(legendreHighPrecision, numberOfPoints)); + } + + /** + * Creates an integrator of the given order, and whose call to the + * {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction) + * integrate} method will perform an integration on the given interval. + * + * @param numberOfPoints Order of the integration rule. + * @param lowerBound Lower bound of the integration interval. + * @param upperBound Upper bound of the integration interval. + * @return a Gauss-Legendre integrator. + * @throws NotStrictlyPositiveException if number of points is not positive + */ + public GaussIntegrator legendreHighPrecision(int numberOfPoints, + double lowerBound, + double upperBound) + throws NotStrictlyPositiveException { + return new GaussIntegrator(transform(getRule(legendreHighPrecision, numberOfPoints), + lowerBound, upperBound)); + } + + /** + * Creates a Gauss-Hermite integrator of the given order. + * The call to the + * {@link SymmetricGaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction) + * integrate} method will perform a weighted integration on the interval + * \([-\infty, +\infty]\): the computed value is the improper integral of + * \(e^{-x^2}f(x)\) + * where \(f(x)\) is the function passed to the + * {@link SymmetricGaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction) + * integrate} method. + * + * @param numberOfPoints Order of the integration rule. + * @return a Gauss-Hermite integrator. + */ + public SymmetricGaussIntegrator hermite(int numberOfPoints) { + return new SymmetricGaussIntegrator(getRule(hermite, numberOfPoints)); + } + + /** + * @param factory Integration rule factory. + * @param numberOfPoints Order of the integration rule. + * @return the integration nodes and weights. + * @throws NotStrictlyPositiveException if number of points is not positive + * @throws DimensionMismatchException if the elements of the rule pair do not + * have the same length. + */ + private static Pair<double[], double[]> getRule(BaseRuleFactory<? extends Number> factory, + int numberOfPoints) + throws NotStrictlyPositiveException, DimensionMismatchException { + return factory.getRule(numberOfPoints); + } + + /** + * Performs a change of variable so that the integration can be performed + * on an arbitrary interval {@code [a, b]}. + * It is assumed that the natural interval is {@code [-1, 1]}. + * + * @param rule Original points and weights. + * @param a Lower bound of the integration interval. + * @param b Lower bound of the integration interval. + * @return the points and weights adapted to the new interval. + */ + private static Pair<double[], double[]> transform(Pair<double[], double[]> rule, + double a, + double b) { + final double[] points = rule.getFirst(); + final double[] weights = rule.getSecond(); + + // Scaling + final double scale = (b - a) / 2; + final double shift = a + scale; + + for (int i = 0; i < points.length; i++) { + points[i] = points[i] * scale + shift; + weights[i] *= scale; + } + + return new Pair<double[], double[]>(points, weights); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/HermiteRuleFactory.java b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/HermiteRuleFactory.java new file mode 100644 index 0000000..abe3fbe --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/HermiteRuleFactory.java @@ -0,0 +1,177 @@ +/* + * 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.math3.analysis.integration.gauss; + +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.util.Pair; +import org.apache.commons.math3.util.FastMath; + +/** + * Factory that creates a + * <a href="http://en.wikipedia.org/wiki/Gauss-Hermite_quadrature"> + * Gauss-type quadrature rule using Hermite polynomials</a> + * of the first kind. + * Such a quadrature rule allows the calculation of improper integrals + * of a function + * <p> + * \(f(x) e^{-x^2}\) + * </p><p> + * Recurrence relation and weights computation follow + * <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun"> + * Abramowitz and Stegun, 1964</a>. + * </p><p> + * The coefficients of the standard Hermite polynomials grow very rapidly. + * In order to avoid overflows, each Hermite polynomial is normalized with + * respect to the underlying scalar product. + * The initial interval for the application of the bisection method is + * based on the roots of the previous Hermite polynomial (interlacing). + * Upper and lower bounds of these roots are provided by </p> + * <blockquote> + * I. Krasikov, + * <em>Nonnegative quadratic forms and bounds on orthogonal polynomials</em>, + * Journal of Approximation theory <b>111</b>, 31-49 + * </blockquote> + * + * @since 3.3 + */ +public class HermiteRuleFactory extends BaseRuleFactory<Double> { + /** π<sup>1/2</sup> */ + private static final double SQRT_PI = 1.77245385090551602729; + /** π<sup>-1/4</sup> */ + private static final double H0 = 7.5112554446494248286e-1; + /** π<sup>-1/4</sup> √2 */ + private static final double H1 = 1.0622519320271969145; + + /** {@inheritDoc} */ + @Override + protected Pair<Double[], Double[]> computeRule(int numberOfPoints) + throws DimensionMismatchException { + + if (numberOfPoints == 1) { + // Break recursion. + return new Pair<Double[], Double[]>(new Double[] { 0d }, + new Double[] { SQRT_PI }); + } + + // Get previous rule. + // If it has not been computed yet it will trigger a recursive call + // to this method. + final int lastNumPoints = numberOfPoints - 1; + final Double[] previousPoints = getRuleInternal(lastNumPoints).getFirst(); + + // Compute next rule. + final Double[] points = new Double[numberOfPoints]; + final Double[] weights = new Double[numberOfPoints]; + + final double sqrtTwoTimesLastNumPoints = FastMath.sqrt(2 * lastNumPoints); + final double sqrtTwoTimesNumPoints = FastMath.sqrt(2 * numberOfPoints); + + // Find i-th root of H[n+1] by bracketing. + final int iMax = numberOfPoints / 2; + for (int i = 0; i < iMax; i++) { + // Lower-bound of the interval. + double a = (i == 0) ? -sqrtTwoTimesLastNumPoints : previousPoints[i - 1].doubleValue(); + // Upper-bound of the interval. + double b = (iMax == 1) ? -0.5 : previousPoints[i].doubleValue(); + + // H[j-1](a) + double hma = H0; + // H[j](a) + double ha = H1 * a; + // H[j-1](b) + double hmb = H0; + // H[j](b) + double hb = H1 * b; + for (int j = 1; j < numberOfPoints; j++) { + // Compute H[j+1](a) and H[j+1](b) + final double jp1 = j + 1; + final double s = FastMath.sqrt(2 / jp1); + final double sm = FastMath.sqrt(j / jp1); + final double hpa = s * a * ha - sm * hma; + final double hpb = s * b * hb - sm * hmb; + hma = ha; + ha = hpa; + hmb = hb; + hb = hpb; + } + + // Now ha = H[n+1](a), and hma = H[n](a) (same holds for b). + // Middle of the interval. + double c = 0.5 * (a + b); + // P[j-1](c) + double hmc = H0; + // P[j](c) + double hc = H1 * c; + boolean done = false; + while (!done) { + done = b - a <= Math.ulp(c); + hmc = H0; + hc = H1 * c; + for (int j = 1; j < numberOfPoints; j++) { + // Compute H[j+1](c) + final double jp1 = j + 1; + final double s = FastMath.sqrt(2 / jp1); + final double sm = FastMath.sqrt(j / jp1); + final double hpc = s * c * hc - sm * hmc; + hmc = hc; + hc = hpc; + } + // Now h = H[n+1](c) and hm = H[n](c). + if (!done) { + if (ha * hc < 0) { + b = c; + hmb = hmc; + hb = hc; + } else { + a = c; + hma = hmc; + ha = hc; + } + c = 0.5 * (a + b); + } + } + final double d = sqrtTwoTimesNumPoints * hmc; + final double w = 2 / (d * d); + + points[i] = c; + weights[i] = w; + + final int idx = lastNumPoints - i; + points[idx] = -c; + weights[idx] = w; + } + + // If "numberOfPoints" is odd, 0 is a root. + // Note: as written, the test for oddness will work for negative + // integers too (although it is not necessary here), preventing + // a FindBugs warning. + if (numberOfPoints % 2 != 0) { + double hm = H0; + for (int j = 1; j < numberOfPoints; j += 2) { + final double jp1 = j + 1; + hm = -FastMath.sqrt(j / jp1) * hm; + } + final double d = sqrtTwoTimesNumPoints * hm; + final double w = 2 / (d * d); + + points[iMax] = 0d; + weights[iMax] = w; + } + + return new Pair<Double[], Double[]>(points, weights); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionRuleFactory.java b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionRuleFactory.java new file mode 100644 index 0000000..e07cb16 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionRuleFactory.java @@ -0,0 +1,215 @@ +/* + * 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.math3.analysis.integration.gauss; + +import java.math.BigDecimal; +import java.math.MathContext; + +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.util.Pair; + +/** + * Factory that creates Gauss-type quadrature rule using Legendre polynomials. + * In this implementation, the lower and upper bounds of the natural interval + * of integration are -1 and 1, respectively. + * The Legendre polynomials are evaluated using the recurrence relation + * presented in <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun"> + * Abramowitz and Stegun, 1964</a>. + * + * @since 3.1 + */ +public class LegendreHighPrecisionRuleFactory extends BaseRuleFactory<BigDecimal> { + /** Settings for enhanced precision computations. */ + private final MathContext mContext; + /** The number {@code 2}. */ + private final BigDecimal two; + /** The number {@code -1}. */ + private final BigDecimal minusOne; + /** The number {@code 0.5}. */ + private final BigDecimal oneHalf; + + /** + * Default precision is {@link MathContext#DECIMAL128 DECIMAL128}. + */ + public LegendreHighPrecisionRuleFactory() { + this(MathContext.DECIMAL128); + } + + /** + * @param mContext Precision setting for computing the quadrature rules. + */ + public LegendreHighPrecisionRuleFactory(MathContext mContext) { + this.mContext = mContext; + two = new BigDecimal("2", mContext); + minusOne = new BigDecimal("-1", mContext); + oneHalf = new BigDecimal("0.5", mContext); + } + + /** {@inheritDoc} */ + @Override + protected Pair<BigDecimal[], BigDecimal[]> computeRule(int numberOfPoints) + throws DimensionMismatchException { + + if (numberOfPoints == 1) { + // Break recursion. + return new Pair<BigDecimal[], BigDecimal[]>(new BigDecimal[] { BigDecimal.ZERO }, + new BigDecimal[] { two }); + } + + // Get previous rule. + // If it has not been computed yet it will trigger a recursive call + // to this method. + final BigDecimal[] previousPoints = getRuleInternal(numberOfPoints - 1).getFirst(); + + // Compute next rule. + final BigDecimal[] points = new BigDecimal[numberOfPoints]; + final BigDecimal[] weights = new BigDecimal[numberOfPoints]; + + // Find i-th root of P[n+1] by bracketing. + final int iMax = numberOfPoints / 2; + for (int i = 0; i < iMax; i++) { + // Lower-bound of the interval. + BigDecimal a = (i == 0) ? minusOne : previousPoints[i - 1]; + // Upper-bound of the interval. + BigDecimal b = (iMax == 1) ? BigDecimal.ONE : previousPoints[i]; + // P[j-1](a) + BigDecimal pma = BigDecimal.ONE; + // P[j](a) + BigDecimal pa = a; + // P[j-1](b) + BigDecimal pmb = BigDecimal.ONE; + // P[j](b) + BigDecimal pb = b; + for (int j = 1; j < numberOfPoints; j++) { + final BigDecimal b_two_j_p_1 = new BigDecimal(2 * j + 1, mContext); + final BigDecimal b_j = new BigDecimal(j, mContext); + final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext); + + // Compute P[j+1](a) + // ppa = ((2 * j + 1) * a * pa - j * pma) / (j + 1); + + BigDecimal tmp1 = a.multiply(b_two_j_p_1, mContext); + tmp1 = pa.multiply(tmp1, mContext); + BigDecimal tmp2 = pma.multiply(b_j, mContext); + // P[j+1](a) + BigDecimal ppa = tmp1.subtract(tmp2, mContext); + ppa = ppa.divide(b_j_p_1, mContext); + + // Compute P[j+1](b) + // ppb = ((2 * j + 1) * b * pb - j * pmb) / (j + 1); + + tmp1 = b.multiply(b_two_j_p_1, mContext); + tmp1 = pb.multiply(tmp1, mContext); + tmp2 = pmb.multiply(b_j, mContext); + // P[j+1](b) + BigDecimal ppb = tmp1.subtract(tmp2, mContext); + ppb = ppb.divide(b_j_p_1, mContext); + + pma = pa; + pa = ppa; + pmb = pb; + pb = ppb; + } + // Now pa = P[n+1](a), and pma = P[n](a). Same holds for b. + // Middle of the interval. + BigDecimal c = a.add(b, mContext).multiply(oneHalf, mContext); + // P[j-1](c) + BigDecimal pmc = BigDecimal.ONE; + // P[j](c) + BigDecimal pc = c; + boolean done = false; + while (!done) { + BigDecimal tmp1 = b.subtract(a, mContext); + BigDecimal tmp2 = c.ulp().multiply(BigDecimal.TEN, mContext); + done = tmp1.compareTo(tmp2) <= 0; + pmc = BigDecimal.ONE; + pc = c; + for (int j = 1; j < numberOfPoints; j++) { + final BigDecimal b_two_j_p_1 = new BigDecimal(2 * j + 1, mContext); + final BigDecimal b_j = new BigDecimal(j, mContext); + final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext); + + // Compute P[j+1](c) + tmp1 = c.multiply(b_two_j_p_1, mContext); + tmp1 = pc.multiply(tmp1, mContext); + tmp2 = pmc.multiply(b_j, mContext); + // P[j+1](c) + BigDecimal ppc = tmp1.subtract(tmp2, mContext); + ppc = ppc.divide(b_j_p_1, mContext); + + pmc = pc; + pc = ppc; + } + // Now pc = P[n+1](c) and pmc = P[n](c). + if (!done) { + if (pa.signum() * pc.signum() <= 0) { + b = c; + pmb = pmc; + pb = pc; + } else { + a = c; + pma = pmc; + pa = pc; + } + c = a.add(b, mContext).multiply(oneHalf, mContext); + } + } + final BigDecimal nP = new BigDecimal(numberOfPoints, mContext); + BigDecimal tmp1 = pmc.subtract(c.multiply(pc, mContext), mContext); + tmp1 = tmp1.multiply(nP); + tmp1 = tmp1.pow(2, mContext); + BigDecimal tmp2 = c.pow(2, mContext); + tmp2 = BigDecimal.ONE.subtract(tmp2, mContext); + tmp2 = tmp2.multiply(two, mContext); + tmp2 = tmp2.divide(tmp1, mContext); + + points[i] = c; + weights[i] = tmp2; + + final int idx = numberOfPoints - i - 1; + points[idx] = c.negate(mContext); + weights[idx] = tmp2; + } + // If "numberOfPoints" is odd, 0 is a root. + // Note: as written, the test for oddness will work for negative + // integers too (although it is not necessary here), preventing + // a FindBugs warning. + if (numberOfPoints % 2 != 0) { + BigDecimal pmc = BigDecimal.ONE; + for (int j = 1; j < numberOfPoints; j += 2) { + final BigDecimal b_j = new BigDecimal(j, mContext); + final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext); + + // pmc = -j * pmc / (j + 1); + pmc = pmc.multiply(b_j, mContext); + pmc = pmc.divide(b_j_p_1, mContext); + pmc = pmc.negate(mContext); + } + + // 2 / pow(numberOfPoints * pmc, 2); + final BigDecimal nP = new BigDecimal(numberOfPoints, mContext); + BigDecimal tmp1 = pmc.multiply(nP, mContext); + tmp1 = tmp1.pow(2, mContext); + BigDecimal tmp2 = two.divide(tmp1, mContext); + + points[iMax] = BigDecimal.ZERO; + weights[iMax] = tmp2; + } + + return new Pair<BigDecimal[], BigDecimal[]>(points, weights); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreRuleFactory.java b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreRuleFactory.java new file mode 100644 index 0000000..cb0bfec --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreRuleFactory.java @@ -0,0 +1,140 @@ +/* + * 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.math3.analysis.integration.gauss; + +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.util.Pair; + +/** + * Factory that creates Gauss-type quadrature rule using Legendre polynomials. + * In this implementation, the lower and upper bounds of the natural interval + * of integration are -1 and 1, respectively. + * The Legendre polynomials are evaluated using the recurrence relation + * presented in <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun"> + * Abramowitz and Stegun, 1964</a>. + * + * @since 3.1 + */ +public class LegendreRuleFactory extends BaseRuleFactory<Double> { + /** {@inheritDoc} */ + @Override + protected Pair<Double[], Double[]> computeRule(int numberOfPoints) + throws DimensionMismatchException { + + if (numberOfPoints == 1) { + // Break recursion. + return new Pair<Double[], Double[]>(new Double[] { 0d }, + new Double[] { 2d }); + } + + // Get previous rule. + // If it has not been computed yet it will trigger a recursive call + // to this method. + final Double[] previousPoints = getRuleInternal(numberOfPoints - 1).getFirst(); + + // Compute next rule. + final Double[] points = new Double[numberOfPoints]; + final Double[] weights = new Double[numberOfPoints]; + + // Find i-th root of P[n+1] by bracketing. + final int iMax = numberOfPoints / 2; + for (int i = 0; i < iMax; i++) { + // Lower-bound of the interval. + double a = (i == 0) ? -1 : previousPoints[i - 1].doubleValue(); + // Upper-bound of the interval. + double b = (iMax == 1) ? 1 : previousPoints[i].doubleValue(); + // P[j-1](a) + double pma = 1; + // P[j](a) + double pa = a; + // P[j-1](b) + double pmb = 1; + // P[j](b) + double pb = b; + for (int j = 1; j < numberOfPoints; j++) { + final int two_j_p_1 = 2 * j + 1; + final int j_p_1 = j + 1; + // P[j+1](a) + final double ppa = (two_j_p_1 * a * pa - j * pma) / j_p_1; + // P[j+1](b) + final double ppb = (two_j_p_1 * b * pb - j * pmb) / j_p_1; + pma = pa; + pa = ppa; + pmb = pb; + pb = ppb; + } + // Now pa = P[n+1](a), and pma = P[n](a) (same holds for b). + // Middle of the interval. + double c = 0.5 * (a + b); + // P[j-1](c) + double pmc = 1; + // P[j](c) + double pc = c; + boolean done = false; + while (!done) { + done = b - a <= Math.ulp(c); + pmc = 1; + pc = c; + for (int j = 1; j < numberOfPoints; j++) { + // P[j+1](c) + final double ppc = ((2 * j + 1) * c * pc - j * pmc) / (j + 1); + pmc = pc; + pc = ppc; + } + // Now pc = P[n+1](c) and pmc = P[n](c). + if (!done) { + if (pa * pc <= 0) { + b = c; + pmb = pmc; + pb = pc; + } else { + a = c; + pma = pmc; + pa = pc; + } + c = 0.5 * (a + b); + } + } + final double d = numberOfPoints * (pmc - c * pc); + final double w = 2 * (1 - c * c) / (d * d); + + points[i] = c; + weights[i] = w; + + final int idx = numberOfPoints - i - 1; + points[idx] = -c; + weights[idx] = w; + } + // If "numberOfPoints" is odd, 0 is a root. + // Note: as written, the test for oddness will work for negative + // integers too (although it is not necessary here), preventing + // a FindBugs warning. + if (numberOfPoints % 2 != 0) { + double pmc = 1; + for (int j = 1; j < numberOfPoints; j += 2) { + pmc = -j * pmc / (j + 1); + } + final double d = numberOfPoints * pmc; + final double w = 2 / (d * d); + + points[iMax] = 0d; + weights[iMax] = w; + } + + return new Pair<Double[], Double[]>(points, weights); + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/SymmetricGaussIntegrator.java b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/SymmetricGaussIntegrator.java new file mode 100644 index 0000000..7fa4884 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/SymmetricGaussIntegrator.java @@ -0,0 +1,103 @@ +/* + * 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.math3.analysis.integration.gauss; + +import org.apache.commons.math3.analysis.UnivariateFunction; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NonMonotonicSequenceException; +import org.apache.commons.math3.util.Pair; + +/** + * This class's implements {@link #integrate(UnivariateFunction) integrate} + * method assuming that the integral is symmetric about 0. + * This allows to reduce numerical errors. + * + * @since 3.3 + */ +public class SymmetricGaussIntegrator extends GaussIntegrator { + /** + * Creates an integrator from the given {@code points} and {@code weights}. + * The integration interval is defined by the first and last value of + * {@code points} which must be sorted in increasing order. + * + * @param points Integration points. + * @param weights Weights of the corresponding integration nodes. + * @throws NonMonotonicSequenceException if the {@code points} are not + * sorted in increasing order. + * @throws DimensionMismatchException if points and weights don't have the same length + */ + public SymmetricGaussIntegrator(double[] points, + double[] weights) + throws NonMonotonicSequenceException, DimensionMismatchException { + super(points, weights); + } + + /** + * Creates an integrator from the given pair of points (first element of + * the pair) and weights (second element of the pair. + * + * @param pointsAndWeights Integration points and corresponding weights. + * @throws NonMonotonicSequenceException if the {@code points} are not + * sorted in increasing order. + * + * @see #SymmetricGaussIntegrator(double[], double[]) + */ + public SymmetricGaussIntegrator(Pair<double[], double[]> pointsAndWeights) + throws NonMonotonicSequenceException { + this(pointsAndWeights.getFirst(), pointsAndWeights.getSecond()); + } + + /** + * {@inheritDoc} + */ + @Override + public double integrate(UnivariateFunction f) { + final int ruleLength = getNumberOfPoints(); + + if (ruleLength == 1) { + return getWeight(0) * f.value(0d); + } + + final int iMax = ruleLength / 2; + double s = 0; + double c = 0; + for (int i = 0; i < iMax; i++) { + final double p = getPoint(i); + final double w = getWeight(i); + + final double f1 = f.value(p); + final double f2 = f.value(-p); + + final double y = w * (f1 + f2) - c; + final double t = s + y; + + c = (t - s) - y; + s = t; + } + + if (ruleLength % 2 != 0) { + final double w = getWeight(iMax); + + final double y = w * f.value(0d) - c; + final double t = s + y; + + s = t; + } + + return s; + } +} diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/package-info.java b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/package-info.java new file mode 100644 index 0000000..066da30 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/package-info.java @@ -0,0 +1,22 @@ +/* + * 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. + */ +/** + * + * Gauss family of quadrature schemes. + * + */ +package org.apache.commons.math3.analysis.integration.gauss; diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/package-info.java b/src/main/java/org/apache/commons/math3/analysis/integration/package-info.java new file mode 100644 index 0000000..f4df3ba --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/integration/package-info.java @@ -0,0 +1,22 @@ +/* + * 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. + */ +/** + * + * Numerical integration (quadrature) algorithms for univariate real functions. + * + */ +package org.apache.commons.math3.analysis.integration; |