\ Adaptive integration using trapezoidal rule \ with Richardson extrapolation \ Integrate a real function from xa to xb \ Forth Scientific Library Algorithm #19 \ V1.1 @(#)adaptint.seq 1.1 08:53:03 11/3/94 \ Usage: USE( fn.name xa xb err )INTEGRAL \ Examples: \ USE( FSQRT % 0 % 1 % 1.E-3 )INTEGRAL FE. 6.666659e-1 ok \ USE( FSQRT % 0 % 2 % 1.E-4 )INTEGRAL FE. 1.885618e0 ok \ : F1 FDUP FSQRT F* ; ok \ USE( F1 % 0 % 1 % 1.E-3 )INTEGRAL FE. 4.000005e-1 ok \ USE( F1 % 0 % 2 % 1.E-4 )INTEGRAL FE. 2.262741e0 ok \ Programmed by J.V. Noble (from "Scientific FORTH" by JVN) \ ANS Standard Program -- version of 10/5/1994 \ This is an ANS Forth program requiring: \ The FLOAT and FLOAT EXT word sets \ Environmental dependencies: \ Assumes independent floating point stack \ Default: 32-bit floating point precision \ Non STANDARD words: \ S>F : S>F S>D D>F ; ( n --) ( F: -- n) \ F=0 ( puts 0 on fpstack) \ } dereference a one-dimensional array. \ as in A{ I } ( base.adr -- base.adr + offset ) \ % (IMMEDIATE word) converts the following text \ to a floating-point literal \ V: define a function vector \ DEFINES (IMMEDIATE) set a vector, as in \ V: DUMMY ; \ : TEST ( xt -- ) DEFINES DUMMY \ DUMMY ; \ 3 5 ' * TEST . 15 ok \ DFVARIABLE define a double-length fp variable \ FINIT to initialize floating point \ function vectoring \ : USE( ' ; \ : V: CREATE ['] NOOP , DOES> @ EXECUTE ; \ : DEFINES ' >BODY STATE @ \ IF POSTPONE LITERAL POSTPONE ! \ ELSE ! THEN ; IMMEDIATE \ (c) Copyright 1994 Julian V. Noble. Permission is granted \ by the author to use this software for any application provided \ the copyright notice is preserved. \ Example Non STANDARD word definitions: 0 CONSTANT REAL*4 1 CONSTANT REAL*8 \ !!! Note: Make sure the sizes in this table are correct for your system !! CREATE LEN_TAB 1 floats , 8 , 8 , 8 , : #CELLS ( n -- #b) LEN_TAB + @ ; \ 1-dim array: A{ I } leaves address of i'th elt of vector A{ : 1ARRAY ( length cells/float --) CREATE DUP , * ALLOT ; : } ( x[0] n -- x[n]) OVER @ * SWAP CELL+ + ; \ =================== no customization is necessary below here =============== \ Data structures 0 S>F FCONSTANT F=0 4 S>F 3 S>F F/ FCONSTANT F=4/3 : DFVARIABLE CREATE 2 FLOATS ALLOT ; 20 REAL*4 #CELLS 1ARRAY X{ \ to make this program double precision 20 REAL*4 #CELLS 1ARRAY E{ \ replace SF@ and SF! with DF@ and DF! 20 REAL*4 #CELLS 1ARRAY F{ \ REAL*4 at left with REAL*8, and 20 REAL*4 #CELLS 1ARRAY I{ \ FVARIABLE with DFVARIABLE 0 VALUE N FVARIABLE OLD.I FVARIABLE FINAL.I \ Begin definitions proper : )INT ( n --) \ trapezoidal rule \ F" ( F(N) + F(N-1) ) * ( X(N) - X(N-1) ) / 2 " X{ OVER } SF@ X{ OVER 1- } SF@ F- F2/ F{ OVER } SF@ F{ OVER 1- } SF@ F+ F* I{ SWAP 1- } SF! ; V: DUMMY \ dummy function name : INITIALIZE ( xt --) ( F: xa xb eps -- integral) DEFINES DUMMY 1 TO N E{ 0 } SF! X{ 1 } SF! X{ 0 } SF! X{ 0 } SF@ DUMMY F{ 0 } SF! \ F" f(0) = dummy( x(0) ) " X{ 1 } SF@ DUMMY F{ 1 } SF! \ F" f(1) = dummy( x(1) ) " 1 )INT F=0 FINAL.I SF! FINIT ; : CHECK.N N 19 > ABORT" Too many subdivisions!" ; : E/2 E{ N 1- } DUP SF@ F2/ SF! ; : }DOWN ( adr n --) OVER @ >R } DUP R@ + R> MOVE ; : MOVE.DOWN E{ N 1- }DOWN X{ N }DOWN F{ N }DOWN ; : X' \ F" X(N) = ( X(N) + X(N-1) ) / 2 " \ F" F(N) = DUMMY(X(N)) " X{ N } SF@ X{ N 1- } SF@ F+ F2/ FDUP X{ N } SF! DUMMY F{ N } SF! ; : N+1 N 1+ TO N ; : N-2 N 2 - TO N ; : SUBDIVIDE CHECK.N E/2 MOVE.DOWN I{ N 1- } SF@ OLD.I SF! X' N )INT N 1+ )INT ; : CONVERGED? ( F: -- I[N]+I'[N-1]-I[N-1]) ( -- f) \ F" I(N) + IP(N-1) - I(N-1) " I{ N } SF@ I{ N 1- } SF@ F+ OLD.I SF@ F- FDUP FABS E{ N 1- } SF@ F2* F< ; : INTERPOLATE ( F: I[N]+I'[N-1]-I[N-1] -- ) \ F" FINAL.I = ( I(N)+I'(N-1) - OLD.I ) * (4/3) + OLD.I + FINAL.I " F=4/3 F* OLD.I SF@ F+ FINAL.I SF@ F+ \ accumulate FINAL.I SF! ; \ store it : )INTEGRAL ( F: A B ERR -- I[A,B]) ( xt --) INITIALIZE BEGIN N 0> WHILE SUBDIVIDE CONVERGED? N+1 IF INTERPOLATE N-2 ELSE FDROP THEN REPEAT FINAL.I SF@ ;