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 */
017 package org.apache.commons.math.analysis.integration;
018
019 import org.apache.commons.math.ConvergenceException;
020 import org.apache.commons.math.FunctionEvaluationException;
021 import org.apache.commons.math.MathRuntimeException;
022 import org.apache.commons.math.MaxIterationsExceededException;
023 import org.apache.commons.math.analysis.UnivariateRealFunction;
024 import org.apache.commons.math.exception.util.LocalizedFormats;
025 import org.apache.commons.math.util.FastMath;
026
027 /**
028 * Implements the <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html">
029 * Legendre-Gauss</a> quadrature formula.
030 * <p>
031 * Legendre-Gauss integrators are efficient integrators that can
032 * accurately integrate functions with few functions evaluations. A
033 * Legendre-Gauss integrator using an n-points quadrature formula can
034 * integrate exactly 2n-1 degree polynomials.
035 * </p>
036 * <p>
037 * These integrators evaluate the function on n carefully chosen
038 * abscissas in each step interval (mapped to the canonical [-1 1] interval).
039 * The evaluation abscissas are not evenly spaced and none of them are
040 * at the interval endpoints. This implies the function integrated can be
041 * undefined at integration interval endpoints.
042 * </p>
043 * <p>
044 * The evaluation abscissas x<sub>i</sub> are the roots of the degree n
045 * Legendre polynomial. The weights a<sub>i</sub> of the quadrature formula
046 * integrals from -1 to +1 ∫ Li<sup>2</sup> where Li (x) =
047 * ∏ (x-x<sub>k</sub>)/(x<sub>i</sub>-x<sub>k</sub>) for k != i.
048 * </p>
049 * <p>
050 * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $
051 * @since 1.2
052 */
053
054 public class LegendreGaussIntegrator extends UnivariateRealIntegratorImpl {
055
056 /** Abscissas for the 2 points method. */
057 private static final double[] ABSCISSAS_2 = {
058 -1.0 / FastMath.sqrt(3.0),
059 1.0 / FastMath.sqrt(3.0)
060 };
061
062 /** Weights for the 2 points method. */
063 private static final double[] WEIGHTS_2 = {
064 1.0,
065 1.0
066 };
067
068 /** Abscissas for the 3 points method. */
069 private static final double[] ABSCISSAS_3 = {
070 -FastMath.sqrt(0.6),
071 0.0,
072 FastMath.sqrt(0.6)
073 };
074
075 /** Weights for the 3 points method. */
076 private static final double[] WEIGHTS_3 = {
077 5.0 / 9.0,
078 8.0 / 9.0,
079 5.0 / 9.0
080 };
081
082 /** Abscissas for the 4 points method. */
083 private static final double[] ABSCISSAS_4 = {
084 -FastMath.sqrt((15.0 + 2.0 * FastMath.sqrt(30.0)) / 35.0),
085 -FastMath.sqrt((15.0 - 2.0 * FastMath.sqrt(30.0)) / 35.0),
086 FastMath.sqrt((15.0 - 2.0 * FastMath.sqrt(30.0)) / 35.0),
087 FastMath.sqrt((15.0 + 2.0 * FastMath.sqrt(30.0)) / 35.0)
088 };
089
090 /** Weights for the 4 points method. */
091 private static final double[] WEIGHTS_4 = {
092 (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0,
093 (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0,
094 (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0,
095 (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0
096 };
097
098 /** Abscissas for the 5 points method. */
099 private static final double[] ABSCISSAS_5 = {
100 -FastMath.sqrt((35.0 + 2.0 * FastMath.sqrt(70.0)) / 63.0),
101 -FastMath.sqrt((35.0 - 2.0 * FastMath.sqrt(70.0)) / 63.0),
102 0.0,
103 FastMath.sqrt((35.0 - 2.0 * FastMath.sqrt(70.0)) / 63.0),
104 FastMath.sqrt((35.0 + 2.0 * FastMath.sqrt(70.0)) / 63.0)
105 };
106
107 /** Weights for the 5 points method. */
108 private static final double[] WEIGHTS_5 = {
109 (322.0 - 13.0 * FastMath.sqrt(70.0)) / 900.0,
110 (322.0 + 13.0 * FastMath.sqrt(70.0)) / 900.0,
111 128.0 / 225.0,
112 (322.0 + 13.0 * FastMath.sqrt(70.0)) / 900.0,
113 (322.0 - 13.0 * FastMath.sqrt(70.0)) / 900.0
114 };
115
116 /** Abscissas for the current method. */
117 private final double[] abscissas;
118
119 /** Weights for the current method. */
120 private final double[] weights;
121
122 /**
123 * Build a Legendre-Gauss integrator.
124 * @param n number of points desired (must be between 2 and 5 inclusive)
125 * @param defaultMaximalIterationCount maximum number of iterations
126 * @exception IllegalArgumentException if the number of points is not
127 * in the supported range
128 */
129 public LegendreGaussIntegrator(final int n, final int defaultMaximalIterationCount)
130 throws IllegalArgumentException {
131 super(defaultMaximalIterationCount);
132 switch(n) {
133 case 2 :
134 abscissas = ABSCISSAS_2;
135 weights = WEIGHTS_2;
136 break;
137 case 3 :
138 abscissas = ABSCISSAS_3;
139 weights = WEIGHTS_3;
140 break;
141 case 4 :
142 abscissas = ABSCISSAS_4;
143 weights = WEIGHTS_4;
144 break;
145 case 5 :
146 abscissas = ABSCISSAS_5;
147 weights = WEIGHTS_5;
148 break;
149 default :
150 throw MathRuntimeException.createIllegalArgumentException(
151 LocalizedFormats.N_POINTS_GAUSS_LEGENDRE_INTEGRATOR_NOT_SUPPORTED,
152 n, 2, 5);
153 }
154
155 }
156
157 /** {@inheritDoc} */
158 @Deprecated
159 public double integrate(final double min, final double max)
160 throws ConvergenceException, FunctionEvaluationException, IllegalArgumentException {
161 return integrate(f, min, max);
162 }
163
164 /** {@inheritDoc} */
165 public double integrate(final UnivariateRealFunction f, final double min, final double max)
166 throws ConvergenceException, FunctionEvaluationException, IllegalArgumentException {
167
168 clearResult();
169 verifyInterval(min, max);
170 verifyIterationCount();
171
172 // compute first estimate with a single step
173 double oldt = stage(f, min, max, 1);
174
175 int n = 2;
176 for (int i = 0; i < maximalIterationCount; ++i) {
177
178 // improve integral with a larger number of steps
179 final double t = stage(f, min, max, n);
180
181 // estimate error
182 final double delta = FastMath.abs(t - oldt);
183 final double limit =
184 FastMath.max(absoluteAccuracy,
185 relativeAccuracy * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5);
186
187 // check convergence
188 if ((i + 1 >= minimalIterationCount) && (delta <= limit)) {
189 setResult(t, i);
190 return result;
191 }
192
193 // prepare next iteration
194 double ratio = FastMath.min(4, FastMath.pow(delta / limit, 0.5 / abscissas.length));
195 n = FastMath.max((int) (ratio * n), n + 1);
196 oldt = t;
197
198 }
199
200 throw new MaxIterationsExceededException(maximalIterationCount);
201
202 }
203
204 /**
205 * Compute the n-th stage integral.
206 * @param f the integrand function
207 * @param min the lower bound for the interval
208 * @param max the upper bound for the interval
209 * @param n number of steps
210 * @return the value of n-th stage integral
211 * @throws FunctionEvaluationException if an error occurs evaluating the
212 * function
213 */
214 private double stage(final UnivariateRealFunction f,
215 final double min, final double max, final int n)
216 throws FunctionEvaluationException {
217
218 // set up the step for the current stage
219 final double step = (max - min) / n;
220 final double halfStep = step / 2.0;
221
222 // integrate over all elementary steps
223 double midPoint = min + halfStep;
224 double sum = 0.0;
225 for (int i = 0; i < n; ++i) {
226 for (int j = 0; j < abscissas.length; ++j) {
227 sum += weights[j] * f.value(midPoint + halfStep * abscissas[j]);
228 }
229 midPoint += step;
230 }
231
232 return halfStep * sum;
233
234 }
235
236 }