steve
/
factor
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Derivatives without dynamics OR locals

db4
Rex Ford 2008-08-12 11:24:00 -04:00
parent 2271aae7f0
commit 359bff5f15
3 changed files with 181 additions and 63 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
extra/math/derivatives/authors.txt
View File

@ -1 +1,2 @@
Reginald Ford
Reginald Ford
Eduardo Cavazos

66
extra/math/derivatives/derivatives-docs.factor
View File

@ -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"
}
}
} ;

175
extra/math/derivatives/derivatives.factor
View File

@ -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 <float-array> ] 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 <float-array> ] 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 ;