On 03-Sep-07 21:45:40, Alberto Monteiro wrote:
> Rory Winston wrote:
>>
>> I am currently (for pedagogical purposes) writing a simple numerical
>> analysis library in R. I have come unstuck when writing a simple
>> Newton-Raphson implementation, that looks like this:
>>
>> f <- function(x) { 2*cos(x)^2 + 3*sin(x) + 0.5 }
>>
>> root <- newton(f, tol=0.0001, N=20, a=1)
>>
>> My issue is calculating the symbolic derivative of f() inside the
>> newton() function.
>>
> If it's pedagogical, maybe returning to basics could help.
>
> What is f'(x)?
>
> It's the limit of (f(x + h) - f(x)) / h when h tends to zero.
>
> So, do it numerically: take a sufficiently small h and
> compute the limit. h must be small enough
> that h^2 f''(x) is much smaller than h f'(x), but big
> enough that f(x+h) is not f(x)
>
> numerical.derivative <- function(f, x, h = 0.0001)
> {
> # test something
> (f(x + h) - f(x)) / h
> }
>
> Ex:
>
> numerical.derivative(cos, pi) = 5e-05 # should be -sin(pi) = 0
> numerical.derivative(sin, pi) = -1 # ok
> numerical.derivative(exp, 0) = 1.00005 # close enough
> numerical.derivative(sqrt, 0) = 100 # should be Inf
>
> Alberto Monteiro
If you want to "go back to basics", it's worth noting that
for a function which has a derivative at x (which excludes
your sqrt(x) at x=0, despite the result you got above,
since only the 1-sided limit as h-->0+ exists):
(f(x+h/2) - f(h-h/2))/h
is generally a far better approximation to f'(x) than is
(f(x+h) - f(x))/h
since the term in h^2 * f''(x) in the expansion is
automatically eliminated! So the accuracy is O(h^2), not
O(h) as with yur definition.
num.deriv<-function(f,x,h=0.001){(f(x + h/2) - f(x-h/2))/h}
num.deriv(cos, pi)
## [1] 0
num.deriv(sin, pi)
[1] -1
but of course num.deriv(sqrt,0) won't work, since it requires
evaluation of sqrt(-h/2).
Hovever, other comparisons with your definition are possible:
numerical.derivative(sin,pi/3)-0.5
##[1] -4.33021e-05
num.deriv(sin,pi/3)-0.5
##[1] -2.083339e-08
Best wishes,
Ted.
--------------------------------------------------------------------
E-Mail: (Ted Harding) <[EMAIL PROTECTED]>
Fax-to-email: +44 (0)870 094 0861
Date: 04-Sep-07 Time: 00:10:42
------------------------------ XFMail ------------------------------
______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
