[R] Derivative of a Function Expression

[(Ted Harding)]

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) <Ted.Harding at manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 04-Sep-07                                       Time: 00:10:42
------------------------------ XFMail ------------------------------