\ Integration by recursion -- an illustration of the concept. \ Uses 3 point Gaussian integration \ version of January 21st, 2003 - 10:51 \ \ Has a Fortran-ish interface for clarity \ \ --------------------------------------------------- \ (c) Copyright 2003 Julian V. Noble. \ \ Permission is granted by the author to \ \ use this software for any application pro- \ \ vided this copyright notice is preserved. \ \ --------------------------------------------------- FALSE [IF] Algorithm: c = (a+b)/2 integral(a,b) = integral(a,c) + integral(c,b) uses 3 point Gauss-Legendre with Richardson extrapolation Environmental dependencies: FLOAT, separate floating point stack Usage: use( fn.name xa xb err )integral Examples: use( fsqrt 0e 1e 1e-4 )integral fs. 6.66666744641E-1 ok use( fsqrt 0e 1e 1e-6 )integral fs. 6.66666670150E-1 ok : f1 FDUP FSQRT F* ; ok use( f1 0e 1e 1e-4 )integral fs. 3.99994557189E-1 ok use( f1 0e 2e 1e-6 )integral fs. 2.26274169632E0 ok [THEN] MARKER -int include ftran201.f \ be sure to review ftrandoc.txt !! \ points and weights for 3 point Gauss-Legendre integration FVARIABLE w0 f$" w0 = 8 / 9 " FVARIABLE w1 f$" w1 = 5 / 9 " FVARIABLE x1 f$" x1 = sqrt( 9 / 15 )" v: fdummy FVARIABLE xa FVARIABLE xb FVARIABLE x0 FVARIABLE dx : scale ( f: xa xb -- ) f" xb =" f" xa =" f" dx = 0.5 * ( xb - xa )" f" x0 = xa + dx" ; VARIABLE Ntimes \ count of function evaluations : gauss3 ( f: xa xb -- gauss_int) \ 3 point Gauss-Legendre scale f" xa = x0 - dx * x1" f" xb = x0 + dx * x1" f" dx * ( w0 * fdummy(x0) + w1 * ( fdummy( xa ) + fdummy( xb ) ) )" 3 Ntimes +! ; FVARIABLE oldI FVARIABLE newI FVARIABLE deltaI FVARIABLE xa FVARIABLE xb FVARIABLE xc FVARIABLE err : juggle ( f: -- xa xc err/2 xc xb err/2 ) f" xa" f" xc" f" 0.5 * err" FOVER FOVER f" xb" FSWAP ; : LessThan F< ( f: x y --) ( -- flag) ; FVARIABLE finI : adaptive ( f: xa xb err -- result) f" err =" f" xb =" f" xa =" f" xc = 0.5 * ( xa + xb )" f" oldI = gauss3( xa, xb )" f" newI = gauss3( xa, xc ) + gauss3( xc, xb )" f" deltaI = newI - oldI" f" LessThan( abs(deltaI), err)" IF f" finI = finI + newI + deltaI / 63" ELSE juggle RECURSE RECURSE THEN ; : )integral ( xt --) ( f: xa xb err -- integral) defines fdummy \ pass function xt to fdummy f" finI = 0" \ initialize integral 0 Ntimes ! \ initialize count adaptive f" finI" Ntimes @ . ." function calls" ;