\ Integration by recursion -- an illustration of the concept \ Uses 3 point Gaussian integration \ \ --------------------------------------------------- \ (c) Copyright 1998 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) 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 -rint : undefined BL WORD FIND NIP 0= ; \ vectoring: for fwd recursion, or using function names as arguments undefined use( [IF] : use( ' \ state-smart ' for syntactic sugar STATE @ IF POSTPONE LITERAL THEN ; IMMEDIATE ' NOOP CONSTANT 'noop : v: CREATE 'noop , DOES> PERFORM ; \ create dummy def'n : 'dfa ' >BODY ; ( -- data field address) : defines 'dfa STATE @ IF POSTPONE LITERAL POSTPONE ! ELSE ! THEN ; IMMEDIATE [THEN] \ end vectoring \ FVALUE undefined FVALUE [IF] : FVALUE CREATE 0e F, DOES> F@ ; : FTO 'dfa STATE @ IF POSTPONE LITERAL POSTPONE F! ELSE F! THEN ; IMMEDIATE [THEN] \ points and weights for 3 point Gauss-Legendre integration 8e 9e F/ FCONSTANT w0 5e 9e F/ FCONSTANT w1 9e 15e F/ FSQRT FCONSTANT x1 x1 FNEGATE FCONSTANT -x1 v: fdummy FVALUE dx FVALUE : scale ( f: xa xb -- ) FOVER F- F2/ ( f: [xb-xa]/2) FSWAP FOVER F+ ( f: [xb-xa]/2 [xa+xb]/2 ) FTO FTO dx ; : )int ( f: xa xb -- gauss_int) scale fdummy w0 F* \ f(0) * w0 x1 dx F* F+ fdummy \ f(x1) -x1 dx F* F+ fdummy F+ w1 F* \ [f(-x1)+f(x1)] * w1 F+ dx F* ; FVALUE oldI FVALUE newI FVALUE Integral FVALUE xa FVALUE xb FVALUE xc FVALUE err : juggle ( f: -- xa xc err/2 xc xb err/2 ) xa xc err F2/ xc FOVER xb FSWAP ; : storeI Integral F+ FTO Integral ; : adaptive ( f: xa xb err -- ) FTO err FTO xb FTO xa xa xb F+ F2/ FTO xc xa xb )int FTO oldI xa xc )int xc xb )int F+ FTO newI newI oldI F- FDUP FABS err F< IF 63e F/ newI F+ storeI ELSE FDROP juggle RECURSE RECURSE THEN ; : )integral ( xt --) ( f: xa xb err -- integral) defines fdummy 0e FTO Integral adaptive Integral ;