diff --git a/extra/math/derivatives/derivatives-docs.factor b/extra/math/derivatives/derivatives-docs.factor index 23847e82f7..70389f18ad 100644 --- a/extra/math/derivatives/derivatives-docs.factor +++ b/extra/math/derivatives/derivatives-docs.factor @@ -3,7 +3,51 @@ USING: help.markup help.syntax ; IN: math.derivatives HELP: derivative ( x function -- m ) -{ $values { "x" "the x-position on the function" } { "function" "a differentiable function" } } -{ $description "Finds the slope of the tangent line at the given x-position on the given function." } ; +{ $values { "x" "a position on the function" } { "function" "a differentiable function" } } +{ $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." + } +} ; -{ derivative-func } related-words +HELP: fast-derivative ( x function -- m ) +{ $values { "x" "a x-position on the function" } { "function" "a differentiable function" } } +{ $description + "Approximates the slope of the tangent line of the provided function " + "by using a secant line with very near points. This implementation is " + "naive and is only provided because it is used in the much more " + "accurate " { $link derivative } " word. Use this word if accuracy " + "is of no importance." +} ; + +HELP: derivative-func ( function -- der ) +{ $values { "func" "a differentiable function" } { "der" "the derivative" } } +{ $description + "Provides the derivative of the function. The implementation simply " + "attaches the " { $link derivative } " word to the end of the function." +} +{ $examples + { $example + "USING: math.derivatives prettyprint ;" + "[ sq ] derivative-func ." + "[ [ sq ] derivative ]" + } +} ; + +ARTICLE: "derivatives" "The Derivative Toolkit" +"A toolkit for computing the derivative of functions." +{ $subsection derivative } +{ $subsection derivative-func } +{ $subsection fast-derivative } ; +ABOUT: "derivatives" diff --git a/extra/math/derivatives/derivatives.factor b/extra/math/derivatives/derivatives.factor index d92066efaf..f77748d0b5 100644 --- a/extra/math/derivatives/derivatives.factor +++ b/extra/math/derivatives/derivatives.factor @@ -1,10 +1,74 @@ -! Copyright © 2008 Reginald Keith Ford II -! Tool for computing the derivative of a function at a point -USING: kernel math math.points math.function-tools ; +! Copyright (c) 2008 Reginald Ford +! 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 ; IN: math.derivatives -: small-amount ( -- n ) 1.0e-14 ; -: some-more ( x -- y ) small-amount + ; -: some-less ( x -- y ) small-amount - ; -: derivative ( x function -- m ) [ [ some-more ] dip eval ] [ [ some-less ] dip eval ] 2bi slope ; -: derivative-func ( function -- function ) [ derivative ] curry ; \ No newline at end of file +! Ridders' method of a derivative, from the chapter +! "Accurate computation of F'(x) and F'(x)F''(x)", +! From "Advances in Engineering Software, Vol. 4, pp. 75-76 +! \ fast-derivative has been factored out for use by children + +: 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 + +: (derivative) ( x function -- m ) + [ [ h get + ] dip eval ] + [ [ h get - ] dip eval ] + 2bi slope ; inline +: fast-derivative ( x function -- m ) + over epsilon sqrt * h set + (derivative) ; 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 +: fast-derivative-func ( function -- function ) [ fast-derivative ] curry ; inline +