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