\ Regula Falsi -- ANS compatible version of September 10th, 1997 \ -- tested with WinForth v. 3.5 on September 10th, 1997 \ Finds roots of real transcendental functions by hybrid \ secant/binary search method \ Environmental dependencies: \ Separate floating point stack \ ANS FLOAT and FLOAT EXT wordsets \ ANS TOOLS EXT wordsets \ Assumes case-insensitivity of definitions \ STANDARD words are UPPERCASE, words defined herein are lowercase \ Expressions are provided as comments, in FORmula TRANslator \ format ( f" expr " ) for the user's convenience. They have \ been tested with the author's program ftran111.f \ (c) Copyright 1997, 1999 Julian V. Noble. Permission is granted \ by the author to use this software for any application provided \ the copyright notice is preserved. MARKER -falsi \ say "-falsi" to evaporate \ conditional compilation BL PARSE undefined DUP PAD C! PAD CHAR+ SWAP CHARS MOVE PAD FIND NIP 0= [IF] : undefined BL WORD FIND NIP 0= ; [THEN] \ Vectoring wordset (conditionally compile if not present) undefined use( [IF] : use( ' STATE @ IF POSTPONE LITERAL THEN ; IMMEDIATE : v: CREATE ['] NOOP , DOES> @ EXECUTE ; : defines ( xt --) ' ( name ) >BODY STATE @ IF POSTPONE LITERAL POSTPONE ! ELSE ! THEN ; IMMEDIATE [THEN] \ Data structures FVARIABLE a \ f(xa) FVARIABLE b \ f(xb) FVARIABLE xa \ lower end of interval FVARIABLE xb \ upper end of interval FVARIABLE epsilon \ precision v: dummy \ vectored function name \ End data structures : x' ( f: -- x') \ secant extrapolation \ f" xa + (xa - xb) * a / (b - A) " ; xa F@ FDUP xb F@ F- ( f: xa xa-xb ) a F@ b F@ FOVER F- F/ F* F+ ; : ( f: -- ) \ binary search extrapolation \ f" (xa + xb) / 2 " ; xa F@ xb F@ F+ F2/ ; : same-sign? ( f: x y --) ( -- f) F* F0> ; : !end ( f: x --) FDUP dummy FDUP ( f: -- x f[x] f[x] ) a F@ same-sign? IF a F! xa F! ELSE b F! xb F! THEN ; : shrink x' !end !end ; \ combine extrapolations : initialize ( xt --) ( f: lower upper precision --) epsilon F! xb F! xa F! \ store parameters defines dummy \ xt -> DUMMY xa F@ dummy a F! \ compute fn at endpts xb F@ dummy b F! a F@ b F@ SAME-SIGN? ABORT" EVEN # OF ROOTS IN INTERVAL!" ; : converged? ( -- f) \ f" ABS( xa - xb ) < EPSILON " ; xa F@ xb F@ F- FABS EPSILON F@ F< ; : )falsi ( xt --) ( f: upper lower precision --) initialize BEGIN shrink converged? UNTIL ; \ Usage example: \ : f1 ( F: x -- [x-e**-x]) FDUP FNEGATE FEXP F- ; \ use( f1 0e0 1e0 1.e-5 )falsi FS. 5.671432E-1 ok