R help - Sep 2007 - Derivative of a Function Expression

Hi

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. I cant seem to get R to do this...I can of course calculate the
derivative by calling D() with an expression object containing the inner
function definition, but I would like to just define the function once and
then compute the derivative of the existing function. I have tried using
deriv() and as.call(), but I am evidently misusing them, as they dont do
what I want. Does anyone know how I can define a function, say foo, which
manipulates one or more arguments, and then refer to that function later in
my code in order to calculate a (partial) derivative?

Thanks
Rory

	[[alternative HTML version deleted]]

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

The Ryacas package can do that (but the function must be one line
and it can't have brace brackets).  The first yacas call below registers f
with
yacas, then we set up a function to act as a template to hold the
derivative and then we set its body calling yacas again to take the
derivative.

library(Ryacas)
f <- function(x)  2*cos(x)^2 + 3*sin(x) +  0.5
yacas(f) # register f with yacas
Df <- f
body(Df) <- yacas(expression(deriv(f(x))))[[1]]
Df

Here is the output:
> library(Ryacas)
> f <- function(x)  2*cos(x)^2 + 3*sin(x) +  0.5
> yacas(f)
[1] "Starting Yacas!"
expression(TRUE)
> Df <- f
> body(Df) <- yacas(expression(deriv(f(x))))[[1]]
> Df
function (x)
2 * (-2 * sin(x) * cos(x)) + 3 * cos(x)

Also see:

demo("Ryacas-Function")

and the other demos, vignette and home page:
   http://ryacas.googlecode.com



On 9/3/07, Rory Winston <rory.winston at gmail.com>
wrote:
> Hi
>
> 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. I cant seem to get R to do this...I can of course calculate the
> derivative by calling D() with an expression object containing the inner
> function definition, but I would like to just define the function once and
> then compute the derivative of the existing function. I have tried using
> deriv() and as.call(), but I am evidently misusing them, as they dont do
> what I want. Does anyone know how I can define a function, say foo, which
> manipulates one or more arguments, and then refer to that function later in
> my code in order to calculate a (partial) derivative?
>
> Thanks
> Rory
>
>        [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at 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.
>