\ 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: use( fname delta x0 )pval ( f: -- integral) \ Example: \ 15 set-precision ok \ use( fexp 1e0 0e0 )pval fs. 2.11450175075134E0 ok \ \ FVARIABLE t \ : f1 ( f: x -- f1) f" t=" f" (t-1)/(1-t^3)" ; \ use( f1 0.25e0 1e0 )pval fs. 1.67823854224122E-1 ok MARKER -pval : undefined BL WORD FIND NIP 0= ; undefined f" [IF] include ftran201.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 f$" xi1 = 0.238619186083197" f$" xi2 = 0.661209386466265" f$" xi3 = 0.932469514203152" f$" wt1 = 0.467913934572691" f$" wt2 = 0.360761573048139" f$" wt3 = 0.171324492379170" 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 " ;