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