\ Adaptive integration using Simpson's rule \ with Richardson extrapolation \ Integrate a real function from xa to xb \ --------------------------------------------------- \ (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. \ \ --------------------------------------------------- \ Usage: use( fn.name xa xb err )integral \ Examples: \ 10 set-precision ok \ : f1 fdup fsqrt f* ; ok \ use( f1 0e0 1e0 1e-6 )integral cr f. 41 function calls \ .4000000556 ok \ use( f1 0e0 1e0 1e-8 )integral cr f. 101 function calls \ .4000000003 ok \ use( fsqrt 0e0 1e0 1e-8 )integral cr f. 189 function calls \ .6666666623 ok \ Programmed by J.V. Noble; version of February 10th, 2003 \ This is an ANS Forth program requiring: \ The FLOAT and FLOAT EXT word sets \ Environmental dependencies: \ Assumes independent floating point stack \ Uses a FORmula TRANslator for clarity MARKER -int needs ftran201.f needs arrays.f FALSE [IF] Non STANDARD words: 1array create a one-dimensional array Ex: 20 1 FLOATS 1array A{ } dereference a one-dimensional array Ex: A{ I } ( base_adr -- base_adr + offset ) Vectoring words included in ftran111.f v: define a function vector Ex: v: dummy defines set a vector Ex: ' * defines dummy 7 2 dummy . 14 ok : test ( xt -- ) defines dummy dummy ; 3 5 ' * test . 15 ok use( get the xt of a word [THEN] WORDLIST CONSTANT adpt \ create separate wordlist adpt SET-CURRENT \ for quadrature def'ns GET-ORDER adpt SWAP 1+ SET-ORDER \ make adpt findable \ Data structures FVARIABLE half f$" half = 0.5" FVARIABLE third f$" third = 1/3" FVARIABLE c43 f$" C43 = 4/3" FVARIABLE fifteenth f$" fifteenth = 1/15" 40 CONSTANT Nmax 1 FLOATS CONSTANT float_len Nmax 2* long float_len 1array x{ Nmax long float_len 1array E{ Nmax 2* long float_len 1array f{ Nmax long float_len 1array I{ 0 VALUE N FVARIABLE oldI FVARIABLE newI FVARIABLE finI FVARIABLE deltaI \ Begin program : )int ( n --) LOCALS| nn | \ Simpson's rule f" I{ nn 2/ 1- } = ( (f{ nn } + f{ nn 2 - }) * third + f{ nn 1- }*c43 ) * ( x{ nn } - x{ nn 1- } ) " ; v: dummy \ dummy function name VARIABLE Ntimes : initialize ( xt --) ( f: xa xb eps -- integral) LOCALS| xt | 2 TO n xt defines dummy E{ 0 } F! \ initial error x{ 2 } F! x{ 0 } F! \ calculate initial points f" x{ 1 } = ( x{ 0 } + x{ 2 } ) * half" f" f{ 0 } = dummy( x{ 0 } ) " \ f_0, f_1 and f_2 f" f{ 1 } = dummy( x{ 1 } ) " f" f{ 2 } = dummy( x{ 2 } ) " 2 )int \ I_0 f" finI = 0" 3 Ntimes ! ; : check.n N [ Nmax 1- ] LITERAL > ABORT" Too many subdivisions!" ; : E/2 f" E{ N 2/ 1- } = half * E{ N 2/ 2 - } " ; \ eps -> eps/2 : move_down f" x{ N } = x{ N 2 - }" f" f{ N } = f{ N 2 - }" f" x{ N 2 - } = x{ N 3 - }" f" f{ N 2 - } = f{ N 3 - }" ; : new_points f" x{ N 1- } = half * ( x{ N } + x{ N 2 - } ) " f" x{ N 3 - } = half * ( x{ N 4 - } + x{ N 2 - } ) " f" f{ N 1- } = dummy( x{ N 1- } ) " f" f{ N 3 - } = dummy( x{ N 3 - } ) " 2 Ntimes +! ; : N+2 N 2 + TO N ; : N-4 N 4 - TO N ; : subdivide N+2 check.n f" oldI = I{ N 2/ 2 - } " E/2 move_down new_points N )int N 2 - )int ; : LessThan F< ; : converged? ( f: --) ( -- f) f" newI = I{ N 2/ 1- } + I{ N 2/ 2 - } " f" deltaI = newI - oldI " f" LessThan( abs(deltaI), E{ N 2/ 2 - })" ; : interpolate! \ use Richardson interpolation and store f" finI = deltaI * fifteenth + newI + finI " ; FORTH-WORDLIST SET-CURRENT \ next definitions to FORTH FALSE VALUE calls-off : )integral ( f: xa xb err -- I[xa,xb]) ( xt --) initialize BEGIN subdivide converged? IF interpolate! N-4 THEN N 0= UNTIL finI F@ calls-off IF Ntimes @ . ." function calls" THEN ; GET-ORDER NIP 1- SET-ORDER \ hide adpt definitions