001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017package org.apache.commons.math4.legacy.analysis.integration; 018 019import org.apache.commons.math4.legacy.analysis.UnivariateFunction; 020import org.apache.commons.math4.legacy.analysis.integration.gauss.GaussIntegrator; 021import org.apache.commons.math4.legacy.analysis.integration.gauss.GaussIntegratorFactory; 022import org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException; 023import org.apache.commons.math4.legacy.exception.TooManyEvaluationsException; 024import org.apache.commons.math4.legacy.exception.util.LocalizedFormats; 025import org.apache.commons.math4.core.jdkmath.JdkMath; 026 027/** 028 * This algorithm divides the integration interval into equally-sized 029 * sub-interval and on each of them performs a 030 * <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html"> 031 * Legendre-Gauss</a> quadrature. 032 * Because of its <em>non-adaptive</em> nature, this algorithm can 033 * converge to a wrong value for the integral (for example, if the 034 * function is significantly different from zero toward the ends of the 035 * integration interval). 036 * In particular, a change of variables aimed at estimating integrals 037 * over infinite intervals as proposed 038 * <a href="http://en.wikipedia.org/w/index.php?title=Numerical_integration#Integrals_over_infinite_intervals"> 039 * here</a> should be avoided when using this class. 040 * 041 * @since 3.1 042 */ 043 044public class IterativeLegendreGaussIntegrator 045 extends BaseAbstractUnivariateIntegrator { 046 /** Factory that computes the points and weights. */ 047 private static final GaussIntegratorFactory FACTORY 048 = new GaussIntegratorFactory(); 049 /** Number of integration points (per interval). */ 050 private final int numberOfPoints; 051 052 /** 053 * Builds an integrator with given accuracies and iterations counts. 054 * 055 * @param n Number of integration points. 056 * @param relativeAccuracy Relative accuracy of the result. 057 * @param absoluteAccuracy Absolute accuracy of the result. 058 * @param minimalIterationCount Minimum number of iterations. 059 * @param maximalIterationCount Maximum number of iterations. 060 * @throws NotStrictlyPositiveException if minimal number of iterations 061 * or number of points are not strictly positive. 062 * @throws org.apache.commons.math4.legacy.exception.NumberIsTooSmallException 063 * if the maximal number of iterations is smaller than or equal to the 064 * minimal number of iterations. 065 */ 066 public IterativeLegendreGaussIntegrator(final int n, 067 final double relativeAccuracy, 068 final double absoluteAccuracy, 069 final int minimalIterationCount, 070 final int maximalIterationCount) { 071 super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount); 072 if (n <= 0) { 073 throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_POINTS, n); 074 } 075 numberOfPoints = n; 076 } 077 078 /** 079 * Builds an integrator with given accuracies. 080 * 081 * @param n Number of integration points. 082 * @param relativeAccuracy Relative accuracy of the result. 083 * @param absoluteAccuracy Absolute accuracy of the result. 084 * @throws NotStrictlyPositiveException if {@code n < 1}. 085 */ 086 public IterativeLegendreGaussIntegrator(final int n, 087 final double relativeAccuracy, 088 final double absoluteAccuracy) { 089 this(n, relativeAccuracy, absoluteAccuracy, 090 DEFAULT_MIN_ITERATIONS_COUNT, DEFAULT_MAX_ITERATIONS_COUNT); 091 } 092 093 /** 094 * Builds an integrator with given iteration counts. 095 * 096 * @param n Number of integration points. 097 * @param minimalIterationCount Minimum number of iterations. 098 * @param maximalIterationCount Maximum number of iterations. 099 * @throws NotStrictlyPositiveException if minimal number of iterations 100 * is not strictly positive. 101 * @throws org.apache.commons.math4.legacy.exception.NumberIsTooSmallException 102 * if the maximal number of iterations is smaller than or equal to the 103 * minimal number of iterations. 104 * @throws NotStrictlyPositiveException if {@code n < 1}. 105 */ 106 public IterativeLegendreGaussIntegrator(final int n, 107 final int minimalIterationCount, 108 final int maximalIterationCount) { 109 this(n, DEFAULT_RELATIVE_ACCURACY, DEFAULT_ABSOLUTE_ACCURACY, 110 minimalIterationCount, maximalIterationCount); 111 } 112 113 /** {@inheritDoc} */ 114 @Override 115 protected double doIntegrate() { 116 // Compute first estimate with a single step. 117 double oldt = stage(1); 118 119 int n = 2; 120 while (true) { 121 // Improve integral with a larger number of steps. 122 final double t = stage(n); 123 124 // Estimate the error. 125 final double delta = JdkMath.abs(t - oldt); 126 final double limit = 127 JdkMath.max(getAbsoluteAccuracy(), 128 getRelativeAccuracy() * (JdkMath.abs(oldt) + JdkMath.abs(t)) * 0.5); 129 130 // check convergence 131 if (iterations.getCount() + 1 >= getMinimalIterationCount() && 132 delta <= limit) { 133 return t; 134 } 135 136 // Prepare next iteration. 137 final double ratio = JdkMath.min(4, JdkMath.pow(delta / limit, 0.5 / numberOfPoints)); 138 n = JdkMath.max((int) (ratio * n), n + 1); 139 oldt = t; 140 iterations.increment(); 141 } 142 } 143 144 /** 145 * Compute the n-th stage integral. 146 * 147 * @param n Number of steps. 148 * @return the value of n-th stage integral. 149 * @throws TooManyEvaluationsException if the maximum number of evaluations 150 * is exceeded. 151 */ 152 private double stage(final int n) 153 throws TooManyEvaluationsException { 154 // Function to be integrated is stored in the base class. 155 final UnivariateFunction f = new UnivariateFunction() { 156 /** {@inheritDoc} */ 157 @Override 158 public double value(double x) { 159 return computeObjectiveValue(x); 160 } 161 }; 162 163 final double min = getMin(); 164 final double max = getMax(); 165 final double step = (max - min) / n; 166 167 double sum = 0; 168 for (int i = 0; i < n; i++) { 169 // Integrate over each sub-interval [a, b]. 170 final double a = min + i * step; 171 final double b = a + step; 172 final GaussIntegrator g = FACTORY.legendreHighPrecision(numberOfPoints, a, b); 173 sum += g.integrate(f); 174 } 175 176 return sum; 177 } 178}