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.math3.analysis.integration;
018
019 import org.apache.commons.math3.exception.MathIllegalArgumentException;
020 import org.apache.commons.math3.exception.MaxCountExceededException;
021 import org.apache.commons.math3.exception.NotStrictlyPositiveException;
022 import org.apache.commons.math3.exception.NumberIsTooSmallException;
023 import org.apache.commons.math3.exception.TooManyEvaluationsException;
024 import org.apache.commons.math3.exception.util.LocalizedFormats;
025 import org.apache.commons.math3.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 function evaluations. A
033 * Legendre-Gauss integrator using an n-points quadrature formula can
034 * integrate 2n-1 degree polynomials exactly.
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 $Id: LegendreGaussIntegrator.java 1364452 2012-07-22 22:30:01Z erans $
051 * @since 1.2
052 * @deprecated As of 3.1 (to be removed in 4.0). Please use
053 * {@link IterativeLegendreGaussIntegrator} instead.
054 */
055 @Deprecated
056 public class LegendreGaussIntegrator extends BaseAbstractUnivariateIntegrator {
057
058 /** Abscissas for the 2 points method. */
059 private static final double[] ABSCISSAS_2 = {
060 -1.0 / FastMath.sqrt(3.0),
061 1.0 / FastMath.sqrt(3.0)
062 };
063
064 /** Weights for the 2 points method. */
065 private static final double[] WEIGHTS_2 = {
066 1.0,
067 1.0
068 };
069
070 /** Abscissas for the 3 points method. */
071 private static final double[] ABSCISSAS_3 = {
072 -FastMath.sqrt(0.6),
073 0.0,
074 FastMath.sqrt(0.6)
075 };
076
077 /** Weights for the 3 points method. */
078 private static final double[] WEIGHTS_3 = {
079 5.0 / 9.0,
080 8.0 / 9.0,
081 5.0 / 9.0
082 };
083
084 /** Abscissas for the 4 points method. */
085 private static final double[] ABSCISSAS_4 = {
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 FastMath.sqrt((15.0 - 2.0 * FastMath.sqrt(30.0)) / 35.0),
089 FastMath.sqrt((15.0 + 2.0 * FastMath.sqrt(30.0)) / 35.0)
090 };
091
092 /** Weights for the 4 points method. */
093 private static final double[] WEIGHTS_4 = {
094 (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0,
095 (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0,
096 (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0,
097 (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0
098 };
099
100 /** Abscissas for the 5 points method. */
101 private static final double[] ABSCISSAS_5 = {
102 -FastMath.sqrt((35.0 + 2.0 * FastMath.sqrt(70.0)) / 63.0),
103 -FastMath.sqrt((35.0 - 2.0 * FastMath.sqrt(70.0)) / 63.0),
104 0.0,
105 FastMath.sqrt((35.0 - 2.0 * FastMath.sqrt(70.0)) / 63.0),
106 FastMath.sqrt((35.0 + 2.0 * FastMath.sqrt(70.0)) / 63.0)
107 };
108
109 /** Weights for the 5 points method. */
110 private static final double[] WEIGHTS_5 = {
111 (322.0 - 13.0 * FastMath.sqrt(70.0)) / 900.0,
112 (322.0 + 13.0 * FastMath.sqrt(70.0)) / 900.0,
113 128.0 / 225.0,
114 (322.0 + 13.0 * FastMath.sqrt(70.0)) / 900.0,
115 (322.0 - 13.0 * FastMath.sqrt(70.0)) / 900.0
116 };
117
118 /** Abscissas for the current method. */
119 private final double[] abscissas;
120
121 /** Weights for the current method. */
122 private final double[] weights;
123
124 /**
125 * Build a Legendre-Gauss integrator with given accuracies and iterations counts.
126 * @param n number of points desired (must be between 2 and 5 inclusive)
127 * @param relativeAccuracy relative accuracy of the result
128 * @param absoluteAccuracy absolute accuracy of the result
129 * @param minimalIterationCount minimum number of iterations
130 * @param maximalIterationCount maximum number of iterations
131 * @exception NotStrictlyPositiveException if minimal number of iterations
132 * is not strictly positive
133 * @exception NumberIsTooSmallException if maximal number of iterations
134 * is lesser than or equal to the minimal number of iterations
135 */
136 public LegendreGaussIntegrator(final int n,
137 final double relativeAccuracy,
138 final double absoluteAccuracy,
139 final int minimalIterationCount,
140 final int maximalIterationCount)
141 throws NotStrictlyPositiveException, NumberIsTooSmallException {
142 super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
143 switch(n) {
144 case 2 :
145 abscissas = ABSCISSAS_2;
146 weights = WEIGHTS_2;
147 break;
148 case 3 :
149 abscissas = ABSCISSAS_3;
150 weights = WEIGHTS_3;
151 break;
152 case 4 :
153 abscissas = ABSCISSAS_4;
154 weights = WEIGHTS_4;
155 break;
156 case 5 :
157 abscissas = ABSCISSAS_5;
158 weights = WEIGHTS_5;
159 break;
160 default :
161 throw new MathIllegalArgumentException(
162 LocalizedFormats.N_POINTS_GAUSS_LEGENDRE_INTEGRATOR_NOT_SUPPORTED,
163 n, 2, 5);
164 }
165
166 }
167
168 /**
169 * Build a Legendre-Gauss integrator with given accuracies.
170 * @param n number of points desired (must be between 2 and 5 inclusive)
171 * @param relativeAccuracy relative accuracy of the result
172 * @param absoluteAccuracy absolute accuracy of the result
173 */
174 public LegendreGaussIntegrator(final int n,
175 final double relativeAccuracy,
176 final double absoluteAccuracy) {
177 this(n, relativeAccuracy, absoluteAccuracy,
178 DEFAULT_MIN_ITERATIONS_COUNT, DEFAULT_MAX_ITERATIONS_COUNT);
179 }
180
181 /**
182 * Build a Legendre-Gauss integrator with given iteration counts.
183 * @param n number of points desired (must be between 2 and 5 inclusive)
184 * @param minimalIterationCount minimum number of iterations
185 * @param maximalIterationCount maximum number of iterations
186 * @exception NotStrictlyPositiveException if minimal number of iterations
187 * is not strictly positive
188 * @exception NumberIsTooSmallException if maximal number of iterations
189 * is lesser than or equal to the minimal number of iterations
190 */
191 public LegendreGaussIntegrator(final int n,
192 final int minimalIterationCount,
193 final int maximalIterationCount) {
194 this(n, DEFAULT_RELATIVE_ACCURACY, DEFAULT_ABSOLUTE_ACCURACY,
195 minimalIterationCount, maximalIterationCount);
196 }
197
198 /** {@inheritDoc} */
199 @Override
200 protected double doIntegrate()
201 throws TooManyEvaluationsException, MaxCountExceededException {
202
203 // compute first estimate with a single step
204 double oldt = stage(1);
205
206 int n = 2;
207 while (true) {
208
209 // improve integral with a larger number of steps
210 final double t = stage(n);
211
212 // estimate error
213 final double delta = FastMath.abs(t - oldt);
214 final double limit =
215 FastMath.max(getAbsoluteAccuracy(),
216 getRelativeAccuracy() * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5);
217
218 // check convergence
219 if ((iterations.getCount() + 1 >= getMinimalIterationCount()) && (delta <= limit)) {
220 return t;
221 }
222
223 // prepare next iteration
224 double ratio = FastMath.min(4, FastMath.pow(delta / limit, 0.5 / abscissas.length));
225 n = FastMath.max((int) (ratio * n), n + 1);
226 oldt = t;
227 iterations.incrementCount();
228
229 }
230
231 }
232
233 /**
234 * Compute the n-th stage integral.
235 * @param n number of steps
236 * @return the value of n-th stage integral
237 * @throws TooManyEvaluationsException if the maximum number of evaluations
238 * is exceeded.
239 */
240 private double stage(final int n)
241 throws TooManyEvaluationsException {
242
243 // set up the step for the current stage
244 final double step = (getMax() - getMin()) / n;
245 final double halfStep = step / 2.0;
246
247 // integrate over all elementary steps
248 double midPoint = getMin() + halfStep;
249 double sum = 0.0;
250 for (int i = 0; i < n; ++i) {
251 for (int j = 0; j < abscissas.length; ++j) {
252 sum += weights[j] * computeObjectiveValue(midPoint + halfStep * abscissas[j]);
253 }
254 midPoint += step;
255 }
256
257 return halfStep * sum;
258
259 }
260
261 }