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 }