\ Cauchy principal value integration using \ 6-point Gauss-Legendre quadrature \ \ --------------------------------------------------- \ (c) Copyright 1999 Julian V. Noble. \ \ Permission is granted by the author to \ \ use this software for any application pro- \ \ vided this copyright notice is preserved. \ \ --------------------------------------------------- \ This is an ANS Forth program requiring the \ FLOAT, FLOAT EXT, FILE and TOOLS EXT wordsets. \ \ Environmental dependencies: \ Assumes independent floating point stack \ Uses a FORmula TRANslator for clarity \ \ Description of algorithm: \ integrate f(x)/(x - x0) from x0-delta to x0+delta \ using even-order Gauss-Legendre quadrature formula \ Usage example: \ 15 set-precision ok \ use( fexp 1e0 0e0 )pval fs. 2.11450175075134E0 ok MARKER -pval : undefined BL WORD FIND NIP 0= ; undefined f" [IF] include ftran111.f [THEN] FVARIABLE wt1 FVARIABLE xi1 FVARIABLE y1 FVARIABLE y1m FVARIABLE wt2 FVARIABLE xi2 FVARIABLE y2 FVARIABLE y2m FVARIABLE wt3 FVARIABLE xi3 FVARIABLE y3 FVARIABLE y3m FVARIABLE x0 FVARIABLE delta 0.238619186083197e0 xi1 F! 0.661209386466265e0 xi2 F! 0.932469514203152e0 xi3 F! 0.467913934572691e0 wt1 F! 0.360761573048139e0 wt2 F! 0.171324492379170e0 wt3 F! v: fdummy : )pval ( xt -- ) ( f: delta x0 -- integral) defines fdummy \ vector function name x0 F! delta F! f" y1 = x0 + delta*xi1" \ scale interval f" y2 = x0 + delta*xi2" f" y3 = x0 + delta*xi3" f" y1m = x0 - delta*xi1" f" y2m = x0 - delta*xi2" f" y3m = x0 - delta*xi3" f" wt1 * ( fdummy(y1) - fdummy(y1m) ) / xi1 + wt2 * ( fdummy(y2) - fdummy(y2m) ) / xi2 + wt3 * ( fdummy(y3) - fdummy(y3m) ) / xi3 " ;