From 359bff5f154caaaaa9ed5037e78d43033cdf48b3 Mon Sep 17 00:00:00 2001 From: Rex Ford Date: Tue, 12 Aug 2008 11:24:00 -0400 Subject: [PATCH] Derivatives without dynamics OR locals --- extra/math/derivatives/authors.txt | 3 +- .../math/derivatives/derivatives-docs.factor | 66 ++++++- extra/math/derivatives/derivatives.factor | 175 ++++++++++++------ 3 files changed, 181 insertions(+), 63 deletions(-) diff --git a/extra/math/derivatives/authors.txt b/extra/math/derivatives/authors.txt index 137b1605da..3be8a6d4d3 100644 --- a/extra/math/derivatives/authors.txt +++ b/extra/math/derivatives/authors.txt @@ -1 +1,2 @@ -Reginald Ford \ No newline at end of file +Reginald Ford +Eduardo Cavazos \ No newline at end of file diff --git a/extra/math/derivatives/derivatives-docs.factor b/extra/math/derivatives/derivatives-docs.factor index 0db52adfa5..a78a697b76 100644 --- a/extra/math/derivatives/derivatives-docs.factor +++ b/extra/math/derivatives/derivatives-docs.factor @@ -1,5 +1,4 @@ -USING: help.markup help.syntax ; - +USING: help.markup help.syntax math.functions ; IN: math.derivatives HELP: derivative ( x function -- m ) @@ -21,6 +20,46 @@ HELP: derivative ( x function -- m ) } } ; +HELP: (derivative) ( x function h err -- m ) +{ $values + { "x" "a position on the function" } + { "function" "a differentiable function" } + { + "h" "distance between the points of the first secant line used for " + "approximation of the tangent. This distance will be divided " + "constantly, by " { $link con } ". See " { $link init-hh } + " for the code which enforces this. H should be .001 to .5 -- too " + "small can cause bad convergence. Also, h should be small enough " + "to give the correct sgn(f'(x)). In other words, if you're expecting " + "a positive derivative, make h small enough to give the same " + "when plugged into the academic limit definition of a derivative. " + "See " { $link update-hh } " for the code which performs this task." + } + { + "err" "maximum tolerance of increase in error. For example, if this " + "is set to 2.0, the program will terminate with its nearest answer " + "when the error multiplies by 2. See " { $link check-safe } " for " + "the enforcing code." + } +} +{ $description + "Approximates the slope of the tangent line by using Ridders' " + "method of computing derivatives, from the chapter \"Accurate computation " + "of F'(x) and F'(x)F''(x)\", from \"Advances in Engineering Software, " + "Vol. 4, pp. 75-76 ." +} +{ $examples + { $example + "USING: math.derivatives prettyprint ;" + "[ sq ] 4 derivative ." + "8" + } + { $notes + "For applied scientists, you may play with the settings " + "in the source file to achieve arbitrary accuracy. " + } +} ; + HELP: derivative-func ( function -- der ) { $values { "func" "a differentiable function" } { "der" "the derivative" } } { $description @@ -30,8 +69,27 @@ HELP: derivative-func ( function -- der ) { $examples { $example "USING: math.derivatives prettyprint ;" - "[ sq ] derivative-func ." - "[ [ sq ] derivative ]" + "60 deg>rad [ sin ] derivative-func call ." + "0.5000000000000173" + } + { $notes + "Without a heavy algebraic system, derivatives must be " + "approximated. With the current settings, there is a fair trade of " + "speed and accuracy; the first 12 digits " + "will always be correct with " { $link sin } " and " { $link cos } + ". The following code performs a minumum and maximum error test." + { $code + "USING: kernel math math.functions math.trig sequences sequences.lib ;" + "360" + "[" + " deg>rad" + " [ [ sin ] derivative-func call ]" + " ! Note: the derivative of sin is cos" + " [ cos ]" + " bi - abs" + "] map minmax" + + } } } ; diff --git a/extra/math/derivatives/derivatives.factor b/extra/math/derivatives/derivatives.factor index 96d0fc3a81..ad8d944bfe 100644 --- a/extra/math/derivatives/derivatives.factor +++ b/extra/math/derivatives/derivatives.factor @@ -1,64 +1,123 @@ -! Tools for approximating derivatives -USING: kernel math math.functions locals generalizations float-arrays sequences -math.constants namespaces math.function-tools math.points math.ranges math.order ; +USING: kernel continuations combinators sequences math + math.order math.ranges accessors float-arrays ; + IN: math.derivatives +TUPLE: state x func h err i j errt fac hh ans a done ; + : largest-float ( -- x ) HEX: 7fefffffffffffff bits>double ; foldable -: ntab 10 ; ! max size of tableau (main accuracy setting) -: con 1.41 ; ! stepsize is decreased by this per-iteration -: con2 1.9881 ; ! con^2 -: initial-h 0.02 ; ! distance of the 2 points of the first secant line -: safe 2.0 ; ! return when current err is SAFE worse than the best - ! \ safe probably should not be changed -SYMBOL: i -SYMBOL: j -SYMBOL: err -SYMBOL: errt -SYMBOL: fac -SYMBOL: h -SYMBOL: ans -SYMBOL: matrix +: ntab ( -- val ) 8 ; +: con ( -- val ) 1.6 ; +: con2 ( -- val ) con con * ; +: big ( -- val ) largest-float ; +: safe ( -- val ) 2.0 ; -: (derivative) ( x function -- m ) - [ [ h get + ] dip eval ] - [ [ h get - ] dip eval ] - 2bi slope ; inline -: init-matrix ( -- ) - ntab [ ntab ] replicate - matrix set ; -: m-set ( value j i -- ) matrix get nth set-nth ; -: m-get ( j i -- n ) matrix get nth nth ; -:: derivative ( x func -- m ) - init-matrix - initial-h h set - x func (derivative) 0 0 m-set - largest-float err set - ntab 1 - [1,b] [| i | - h [ con / ] change - x func (derivative) 0 i m-set - con2 fac set - i [1,b] [| j | - j 1 - i m-get fac get * - j 1 - i 1 - m-get - - - fac get 1 - - / j i m-set - fac [ con2 * ] change - j i m-get j 1 - i m-get - abs - j i m-get j 1 - i 1 - m-get - abs - max errt set - errt get err get <= - [ - errt get err set - j i m-get ans set - ] [ ] - if - ] each - i i m-get i 1 - dup m-get - abs - err get safe * - < - ] all? drop - ans get ; inline -: derivative-func ( function -- function ) [ derivative ] curry ; inline +! Yes, this was ported from C code. +: a[i][i] ( state -- elt ) [ i>> ] [ i>> ] [ a>> ] tri nth nth ; +: a[j][i] ( state -- elt ) [ i>> ] [ j>> ] [ a>> ] tri nth nth ; +: a[j-1][i] ( state -- elt ) [ i>> ] [ j>> 1 - ] [ a>> ] tri nth nth ; +: a[j-1][i-1] ( state -- elt ) [ i>> 1 - ] [ j>> 1 - ] [ a>> ] tri nth nth ; +: a[i-1][i-1] ( state -- elt ) [ i>> 1 - ] [ i>> 1 - ] [ a>> ] tri nth nth ; +: check-h ( state -- state ) + dup h>> 0 = [ "h must be nonzero in dfridr" throw ] when ; +: init-a ( state -- state ) ntab [ ntab ] replicate >>a ; +: init-hh ( state -- state ) dup h>> >>hh ; +: init-err ( state -- state ) big >>err ; +: update-hh ( state -- state ) dup hh>> con / >>hh ; +: reset-fac ( state -- state ) con2 >>fac ; +: update-fac ( state -- state ) dup fac>> con2 * >>fac ; + +! If error is decreased, save the improved answer +: error-decreased? ( state -- state ? ) [ ] [ errt>> ] [ err>> ] tri <= ; +: save-improved-answer ( state -- state ) + dup err>> >>errt + dup a[j][i] >>ans ; + +! If higher order is worse by a significant factor SAFE, then quit early. +: check-safe ( state -- state ) + dup + [ [ a[i][i] ] [ a[i-1][i-1] ] bi - abs ] [ err>> safe * ] bi >= + [ t >>done ] + when ; +: x+hh ( state -- val ) [ x>> ] [ hh>> ] bi + ; +: x-hh ( state -- val ) [ x>> ] [ hh>> ] bi - ; +: limit-approx ( state -- val ) + [ + [ [ x+hh ] [ func>> ] bi call ] + [ [ x-hh ] [ func>> ] bi call ] + bi - + ] + [ hh>> 2.0 * ] + bi / ; +: a[0][0]! ( state -- state ) + { [ ] [ limit-approx ] [ drop 0 ] [ drop 0 ] [ a>> ] } cleave nth set-nth ; +: a[0][i]! ( state -- state ) + { [ ] [ limit-approx ] [ i>> ] [ drop 0 ] [ a>> ] } cleave nth set-nth ; +: a[j-1][i]*fac ( state -- val ) [ a[j-1][i] ] [ fac>> ] bi * ; +: new-a[j][i] ( state -- val ) + [ [ a[j-1][i]*fac ] [ a[j-1][i-1] ] bi - ] + [ fac>> 1.0 - ] + bi / ; +: a[j][i]! ( state -- state ) + { [ ] [ new-a[j][i] ] [ i>> ] [ j>> ] [ a>> ] } cleave nth set-nth ; + +: update-errt ( state -- state ) + dup + [ [ a[j][i] ] [ a[j-1][i] ] bi - abs ] + [ [ a[j][i] ] [ a[j-1][i-1] ] bi - abs ] + bi max + >>errt ; + +: not-done? ( state -- state ? ) dup done>> not ; + +: derive ( state -- state ) + init-a + check-h + init-hh + a[0][0]! + init-err + 1 ntab [a,b) + [ + >>i + not-done? + [ + update-hh + a[0][i]! + reset-fac + 1 over i>> [a,b] + [ + >>j + a[j][i]! + update-fac + update-errt + error-decreased? [ save-improved-answer ] when + ] + each + check-safe + ] + when + ] + each ; + +: derivative-state ( x func h err -- state ) + state new + swap >>err + swap >>h + swap >>func + swap >>x ; + +! For scientists: +! h should be .001 to .5 -- too small can cause bad convergence, +! h should be small enough to give the correct sgn(f'(x)) +! err is the max tolerance of gain in error for a single iteration- +: (derivative) ( x func h err -- ans error ) + derivative-state + derive + [ ans>> ] + [ errt>> ] + bi ; + +: derivative ( x func -- m ) 0.01 2.0 (derivative) drop ; +: derivative-func ( func -- der ) [ derivative ] curry ; \ No newline at end of file