001/* 002* Quadrature -- Methods for integrating a 1D function: I=int_a^b{f(x) dx}. 003* 004* Copyright (C) 2009-2014 by Joseph A. Huwaldt. 005* All rights reserved. 006* 007* This library is free software; you can redistribute it and/or 008* modify it under the terms of the GNU Lesser General Public 009* License as published by the Free Software Foundation; either 010* version 2.1 of the License, or (at your option) any later version. 011* 012* This library is distributed in the hope that it will be useful, 013* but WITHOUT ANY WARRANTY; without even the implied warranty of 014* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 015* Lesser General Public License for more details. 016* 017* You should have received a copy of the GNU Lesser General Public License 018* along with this program; if not, write to the Free Software 019* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 020* Or visit: http://www.gnu.org/licenses/lgpl.html 021*/ 022package jahuwaldt.tools.math; 023 024import static java.lang.Math.*; 025 026 027/** 028* A collection of methods for integrating a 1D function within a range: 029* <code>I=int_a^b{f(x) dx}</code>. 030* 031* <p> Modified by: Joseph A. Huwaldt </p> 032* 033* @author Joseph A. Huwaldt Date: June 8, 2009 034* @version January 23, 2014 035*/ 036public class Quadrature { 037 038 // 2^(JMAX-1) is the maximum number of steps in simpsonsRule. 039 private static final int JMAX = 20; 040 041 /** 042 * Returns the integral of the supplied function from <code>a</code> to 043 * <code>b</code> using Simpson's Rule. Reference: Numerical Recipes in C, pg. 139. 044 * 045 * @param func The function being integrated (the "derivative" method is never called). 046 * @param a The lower integration limit (must be < b). 047 * @param b The upper integration limit (must be > a). 048 * @param tol The desired relative error. 049 * @return The integral of the supplied function from <code>a</code> to <code>b</code>. 050 * @throws IntegratorException if to many iterations are attempted. 051 * @throws RootException thrown by the function itself. 052 */ 053 public static double simpsonsRule(Evaluatable1D func, double a, double b, double tol) throws IntegratorException, RootException { 054 if (tol < MathTools.EPS) tol = MathTools.EPS; 055 056 double ost = 0; 057 double os = 0; 058 for (int j = 1; j <= JMAX; ++j) { 059 double st = trapzd(func, a, b, j, ost); 060 double s = (4 * st - ost) / 3.; 061 if (j > 5) 062 if (abs(s - os) < tol * abs(os) || (s == 0.0 && os == 0.0)) 063 return s; 064 os = s; 065 ost = st; 066 } 067 throw new IntegratorException("To many steps in simpsonsRule."); 068 } 069 070 /** 071 * This routine computes the nth stage of refinement of an extended 072 * trapezoidal rule. When called with n=1, the routine returns the 073 * crudest estimate of int_a^b{f(x)dx}. Subsequent calls with 074 * n=2,3,... (in sequential order) will improve the accuracy by adding 075 * 2^(n-2) additional interior points. 076 * Reference: Numerical Recipes in C, pg. 137. 077 * 078 * @param func The function being integrated (the "derivative" function is never called). 079 * @param a The lower integration limit (must be < b). 080 * @param b The upper integration limit (must be > a). 081 * @param n Desired stage of refinement (n=1, 2, 3, ... in sequential order). 082 * @param s The estimate of the integral from the last call to this method (ignored when n==1). 083 * @return The estimate of the function's integral. 084 * @throws RootException thrown by the function itself. 085 */ 086 private static double trapzd(Evaluatable1D func, double a, double b, int n, double s) throws RootException { 087 if (n == 1) return 0.5*(b-a)*(func.function(a)+func.function(b)); 088 089 int it = 1; 090 for (int j=1; j < n-1; j++) it <<= 1; 091 092 double tnm = it; 093 double del = (b-a)/tnm; 094 double x = a + 0.5*del; 095 double sum = 0.; 096 for (int j=0; j < it; j++, x += del) 097 sum += func.function(x); 098 099 s = 0.5*(s+(b-a)*sum/tnm); 100 101 return s; 102 } 103 104 105 /** 106 * Returns the integral of the supplied function from <code>a</code> to <code>b</code> 107 * using Gauss-Legendre integration (W(x) = 1). 108 * Reference: Numerical Recipes in C, pg. 148. 109 * 110 * @param func The function being integrated (the "derivative" function is never called). 111 * @param a The lower integration limit (must be < b). 112 * @param b The upper integration limit (must be > a). 113 * @param x The abscissas for the quadrature. Number of points determined by length of x. 114 * @param w The weights for the quadrature. Must be same length as x. 115 * @return The integral of the supplied function from <code>a</code> to <code>b</code>. 116 * @throws RootException may be thrown by the function itself. 117 */ 118 private static double gaussLegendre_Wx1(Evaluatable1D func, double a, double b, double[] x, double[] w) throws RootException { 119 120 double xm = 0.5*(b+a); 121 double xr = 0.5*(b-a); 122 double s = 0; 123 int size = x.length; 124 for (int j=0; j < size; j++) { 125 double dx = xr*x[j]; 126 s += w[j]*(func.function(xm+dx) + func.function(xm-dx)); 127 } 128 return s *= xr; 129 } 130 131 // Abscissa used by gaussLegendre_Wx1N10(). 132 private static final double[] _gaussLx_Wx1N10 = 133 {0.14887433898163122, 0.43339539412924716, 0.6794095682990244, 0.8650633666889845, 0.9739065285171717}; 134 // Weights used by gaussLegendre_Wx1N10(). 135 private static final double[] _gaussLw_Wx1N10 = 136 {0.2955242247147529, 0.26926671930999174, 0.21908636251598218, 0.14945134915058053, 0.06667134430868371}; 137 138 /** 139 * Returns the integral of the supplied function from <code>a</code> to <code>b</code> 140 * using Gauss-Legendre integration (W(x) = 1, N=10). 141 * Reference: Numerical Recipes in C, pg. 148. 142 * 143 * @param func The function being integrated (the "derivative" function is never called). 144 * @param a The lower integration limit (must be < b). 145 * @param b The upper integration limit (must be > a). 146 * @return The integral of the supplied function from <code>a</code> to <code>b</code>. 147 * @throws RootException may be thrown by the function itself. 148 */ 149 public static double gaussLegendre_Wx1N10(Evaluatable1D func, double a, double b) throws RootException { 150 return gaussLegendre_Wx1(func, a, b, _gaussLx_Wx1N10, _gaussLw_Wx1N10); 151 } 152 153 // Abscissa used by gaussLegendre_Wx1N20(). 154 private static final double[] _gaussLx_Wx1N20 = 155 {0.07652652113349734, 0.2277858511416451, 0.37370608871541955, 0.5108670019508271, 0.636053680726515, 156 0.7463319064601508, 0.8391169718222189, 0.912234428251326, 0.9639719272779138, 0.9931285991850949}; 157 // Weights used by gaussLegendre_Wx1N20(). 158 private static final double[] _gaussLw_Wx1N20 = 159 {0.15275338713072598, 0.14917298647260382, 0.1420961093183819, 0.1316886384491766, 0.11819453196151775, 160 0.10193011981724048, 0.08327674157670474, 0.06267204833410904, 0.04060142980038705, 0.017614007139150577}; 161 162 /** 163 * Returns the integral of the supplied function from <code>a</code> to <code>b</code> 164 * using Gauss-Legendre integration (W(x) = 1, N=20). 165 * Reference: Numerical Recipes in C, pg. 148. 166 * 167 * @param func The function being integrated (the "derivative" function is never called). 168 * @param a The lower integration limit (must be < b). 169 * @param b The upper integration limit (must be > a). 170 * @return The integral of the supplied function from <code>a</code> to <code>b</code>. 171 * @throws RootException may be thrown by the function itself. 172 */ 173 public static double gaussLegendre_Wx1N20(Evaluatable1D func, double a, double b) throws RootException { 174 return gaussLegendre_Wx1(func, a, b, _gaussLx_Wx1N20, _gaussLw_Wx1N20); 175 } 176 177 // Abscissa used by gaussLegendre_Wx1N40(). 178 private static final double[] _gaussLx_Wx1N40 = 179 {0.038772417506050816, 0.11608407067525521, 0.1926975807013711, 0.2681521850072, 0.3419940908257585, 180 0.413779204371605, 0.4830758016861787, 0.5494671250951282, 0.6125538896679802, 0.6719566846141796, 181 0.7273182551899271, 0.7783056514265194, 0.8246122308333117, 0.8659595032122595, 0.9020988069688743, 182 0.9328128082786765, 0.9579168192137917, 0.9772599499837743, 0.990726238699457, 0.9982377097105593}; 183 // Weights used by gaussLegendre_Wx1N40(). 184 private static final double[] _gaussLw_Wx1N40 = 185 {0.0775059479784248, 0.077039818164248, 0.07611036190062617, 0.07472316905796833, 0.07288658239580408, 186 0.07061164739128681, 0.06791204581523393, 0.06480401345660108, 0.06130624249292889, 0.05743976909939157, 187 0.05322784698393679, 0.04869580763507221, 0.04387090818567314, 0.03878216797447161, 0.033460195282546436, 188 0.02793700698001634, 0.02224584919416689, 0.016421058381907876, 0.010498284531152905, 0.004521277098533099}; 189 190 /** 191 * Returns the integral of the supplied function from <code>a</code> to <code>b</code> 192 * using Gauss-Legendre integration (W(x) = 1, N=40). 193 * Reference: Numerical Recipes in C, pg. 148. 194 * 195 * @param func The function being integrated (the "derivative" function is never called). 196 * @param a The lower integration limit (must be < b). 197 * @param b The upper integration limit (must be > a). 198 * @return The integral of the supplied function from <code>a</code> to <code>b</code>. 199 * @throws RootException may be thrown by the function itself. 200 */ 201 public static double gaussLegendre_Wx1N40(Evaluatable1D func, double a, double b) throws RootException { 202 return gaussLegendre_Wx1(func, a, b, _gaussLx_Wx1N40, _gaussLw_Wx1N40); 203 } 204 205 206 // Constants used by the adaptive Lobatto integration routines. 207 private static final double ALPHA = sqrt(2./3.); 208 private static final double BETA = 1./sqrt(5.); 209 private static final double X1 = 0.942882415695480; 210 private static final double X2 = 0.641853342345781; 211 private static final double X3 = 0.236383199662150; 212 213 /** 214 * Returns the integral of the supplied function from <code>a</code> to <code>b</code> 215 * using the adaptive Gauss-Lobatto method.<br> 216 * Reference: 217 * Gander, W., Gautschi, W., "Adaptive Quadrature - Revisited", ETH Zürich, Aug. 1998. 218 * 219 * @param func The function being integrated (the "derivative" function is never called). 220 * @param a The lower integration limit (must be < b). 221 * @param b The upper integration limit (must be > a). 222 * @param tol The desired relative error. 223 * @return The integral of the supplied function from <code>a</code> to <code>b</code>. 224 * @throws IntegratorException if to many iterations are attempted. 225 * @throws RootException may be thrown by the function itself. 226 */ 227 public static double adaptLobatto(Evaluatable1D func, double a, double b, double tol) throws IntegratorException, RootException { 228 if (tol < MathTools.EPS) tol = MathTools.EPS; 229 230 double m = (a+b)/2; 231 double h = (b-a)/2; 232 double[] x = {a, m-X1*h, m-ALPHA*h, m-X2*h, m-BETA*h, m-X3*h, m, m+X3*h, 233 m+BETA*h, m+X2*h, m+ALPHA*h, m+X1*h, b}; 234 double[] y = new double[13]; 235 for (int i=0; i < 13; ++i) 236 y[i] = func.function(x[i]); 237 double fa = y[0]; 238 double fb = y[12]; 239 double i2 = h/6*(fa + fb + 5*(y[4] + y[8])); 240 double i1 = h/1470*(77*(fa + fb) + 432*(y[2] + y[10]) + 625*(y[4] + y[8]) + 672*y[6]); 241 double is = h*(0.01582719197343802*(fa+fb) + 0.0942738402188500*(y[1]+y[11]) 242 + 0.155071987336585*(y[2]+y[10]) + 0.188821573960182*(y[3]+y[9]) 243 + 0.199773405226859*(y[4]+y[8]) + 0.224926465333340*(y[5]+y[7]) 244 + 0.242611071901408*y[6]); 245 double s = signum(is); 246 if (s == 0) s = 1; 247 double erri1 = abs(i1-is); 248 double erri2 = abs(i2-is); 249 double R = erri1/erri2; 250 if (R > 0 && R < 1) 251 tol = tol/R; 252 is = s*abs(is)*tol/MathTools.EPS; 253 if (is == 0) is = b-a; 254 255 return adaptlobstp(func,a,b,fa,fb,is); 256 } 257 258 /** 259 * Recursion function used by adaptlob(). 260 */ 261 private static double adaptlobstp(Evaluatable1D func, double a, double b, double fa, double fb, double is) throws IntegratorException, RootException { 262 double h = (b-a)/2; 263 double m = (a+b)/2; 264 double mll = m-ALPHA*h; 265 double ml = m-BETA*h; 266 double mr = m+BETA*h; 267 double mrr = m+ALPHA*h; 268 double fmll = func.function(mll); 269 double fml = func.function(ml); 270 double fm = func.function(m); 271 double fmr = func.function(mr); 272 double fmrr = func.function(mrr); 273 double i2 = h/6*(fa+fb+5*(fml+fmr)); 274 double i1 = h/1470*(77*(fa + fb) + 432*(fmll + fmrr) + 625*(fml + fmr) + 672*fm); 275 if (is + (i1-i2) == is || mll <= a || b <= mr) { 276 if (m <= a || b <= m ) { 277 throw new IntegratorException("Interval contains no more machine precision. Required tolerance not met."); 278 } 279 return i1; 280 281 } else { 282 double output = adaptlobstp(func,a,mll,fa,fmll,is); 283 output += adaptlobstp(func,mll,ml,fmll,fml,is); 284 output += adaptlobstp(func,ml,m,fml,fm,is); 285 output += adaptlobstp(func,m,mr,fm,fmr,is); 286 output += adaptlobstp(func,mr,mrr,fmr,fmrr,is); 287 output += adaptlobstp(func,mrr,b,fmrr,fb,is); 288 289 return output; 290 } 291 } 292 293} 294 295 296