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 }