I've been trying to write a program in R that implements Newton's method. I've been mostly successful, but there are two little snags that have been bothering me. Here's my code: Newton<-function(f,f.,guess){ #f <- readline(prompt="Function? ") #f. <- readline(prompt="Derivative? ") #guess <- as.numeric(readline(prompt="Guess? ")) a <- rep(NA, length=1000) a[1] <- guess a[2] <- a[1] - f(a[1]) / f.(a[1]) for(i in 2:length(a)){ if(a[i] == a[i-1]){ break } else{ a[i+1] <- a[i] - f(a[i]) / f.(a[i]) } } a <- a[complete.cases(a)] return(a) } I can't get R to recognize the functions f and f. if I try using readline() to prompt for user input. I get the error "Error in Newton() : could not find function "f."" However, if I comment out the readlines (as above), define f and f. beforehand, then everything works fine. I've been trying to make R calculate the derivative of a function. The problem is that the class object with which R can take symbolic derivatives is expression() , but I want to take the derivative of a function() and have it give me a function() . In short, I'm having trouble with type conversion between expression() and function() . I have an ugly but effective solution for going from function() to expression() . Given a function f, D(body(f)[[2]],"x") will give the derivative of f . However, this output is an expression() , and I haven't been able to turn it back into a function() . Do I need to use eval() or something? I've tried subsetting, but to no avail. For instance: g <- expression(sin(x)) g[[1]] sin(x) f <- function(x){g[[1]]} f(0) sin(x) when what I want is f(0) = 0 since sin(0) = 0. EDIT: Thanks all! Here's my new code: Newton<-function(f,f.,guess){ g<-readline(prompt="Function? ") g<-parse(text=g) g.<-D(g,"x") f<-function(x){eval(g[[1]])} f.<-function(x){eval(g.)} guess<-as.numeric(readline(prompt="Guess? ")) a<-rep(NA, length=1000) a[1]<-guess a[2]<-a[1]-f(a[1])/f.(a[1]) for(i in 2:length(a)){ if(a[i]==a[i-1]){break }else{ a[i+1]<-a[i]-f(a[i])/f.(a[i]) } } a<-a[complete.cases(a)] #a<-a[1:(length(a)-1)] return(a) }

This first problem arises because readline reads in a text string, whereas what you need is an expression. You can use parse() to convert the text string to an expression: f <-readline(prompt="Function? ") sin(x) f # [1] "sin(x)" f <- parse(text = f) f # expression(sin(x)) g <- D(f, "x") g # cos(x) To pass in values for the arguments in the function call in the expression (whew!), you can eval() it in an environment containing the supplied values. Nicely, R will allow you to supply those values in a list supplied to the envir= argument of eval() : > eval(f, envir=list(x=0)) # [1] 0

rfunctionexpressiontype-conversion

Viktor Vaughn

Reputation: 223

R type conversion expression() function()

I've been trying to write a program in R that implements Newton's method. I've been mostly successful, but there are two little snags that have been bothering me. Here's my code:

Newton<-function(f,f.,guess){
    #f <- readline(prompt="Function? ")
    #f. <- readline(prompt="Derivative? ")
    #guess <- as.numeric(readline(prompt="Guess? "))
    a <- rep(NA, length=1000)
    a[1] <- guess
    a[2] <- a[1] - f(a[1]) / f.(a[1])
    for(i in 2:length(a)){
        if(a[i] == a[i-1]){
           break
        } 
        else{
           a[i+1] <- a[i] - f(a[i]) / f.(a[i])
        }
    }   
    a <- a[complete.cases(a)]
    return(a)
}

I can't get R to recognize the functions f and f. if I try using readline() to prompt for user input. I get the error "Error in Newton() : could not find function "f."" However, if I comment out the readlines (as above), define f and f. beforehand, then everything works fine.
I've been trying to make R calculate the derivative of a function. The problem is that the class object with which R can take symbolic derivatives is expression(), but I want to take the derivative of a function() and have it give me a function(). In short, I'm having trouble with type conversion between expression() and function().

I have an ugly but effective solution for going from function() to expression(). Given a function f, D(body(f)[[2]],"x") will give the derivative of f. However, this output is an expression(), and I haven't been able to turn it back into a function(). Do I need to use eval() or something? I've tried subsetting, but to no avail. For instance:

g <- expression(sin(x))
g[[1]]
sin(x)
f <- function(x){g[[1]]}
f(0)
sin(x)

when what I want is f(0) = 0 since sin(0) = 0.

EDIT: Thanks all! Here's my new code:

Newton<-function(f,f.,guess){
    g<-readline(prompt="Function? ")
    g<-parse(text=g)
    g.<-D(g,"x")
    f<-function(x){eval(g[[1]])}
    f.<-function(x){eval(g.)}
    guess<-as.numeric(readline(prompt="Guess? "))
    a<-rep(NA, length=1000)
    a[1]<-guess
    a[2]<-a[1]-f(a[1])/f.(a[1])
    for(i in 2:length(a)){
        if(a[i]==a[i-1]){break
        }else{
        a[i+1]<-a[i]-f(a[i])/f.(a[i])
        }
    }   
a<-a[complete.cases(a)]
#a<-a[1:(length(a)-1)]
return(a)
}

Upvotes: 12

Answers (3)

Carl Witthoft

Reputation: 21532

BTW, having recently written a toy which calculates fractal patterns based on root convergence of Newton's method in the complex plane, I can recommend you toss in some code like the following (where the main function's argument list includes "func" and "varname" ).

func<- gsub(varname, 'zvar', func)
    funcderiv<- try( D(parse(text=func), 'zvar') )
    if(class(funcderiv) == 'try-error') stop("Can't calculate derivative")

If you're more cautious, you could include a an argument "funcderiv" , and wrap my code in

if(missing(funcderiv)){blah blah}

Ahh, why not: here's my complete function for all to use and enjoy:-)

# build Newton-Raphson fractal
#define: f(z)  the convergence per Newton's method is 
# zn+1 = zn - f(zn)/f'(zn)
#record which root each starting z0 converges to, 
# and to get even nicer coloring, record the number of iterations to get there.
# Inputs:
#   func: character string, including the variable. E.g., 'x+ 2*x^2' or 'sin(x)'
#   varname: character string indicating the variable name
#   zreal: vector(preferably) of Re(z)
#   zim: vector of Im(z)
#   rootprec: convergence precision for the NewtonRaphson algorithm
#   maxiter: safety switch, maximum iterations, after which throw an error
#
nrfrac<-function(func='z^5 - 1 ', varname = 'z', zreal= seq(-5,5,by=.1), zim, rootprec=1.0e-5, maxiter=1e4, drawplot=T, drawiterplot=F, ...) {
    zreal<-as.vector(zreal)
    if(missing(zim)) zim <- as.vector(zreal)
# precalculate F/F' 
    # check for differentiability (in R's capability)
    # and make sure to get the correct variable name into the function
    func<- gsub(varname, 'zvar', func)
    funcderiv<- try( D(parse(text=func), 'zvar') )
    if(class(funcderiv) == 'try-error') stop("Can't calculate derivative")  
# Interesting "feature" of deparse : default is to limit each string to 60 or64
# chars.  Need to avoid that here.  Doubt I'd ever see a derivative w/ more
# than 500 chars, the max allowed by deparse. To do it right, 
# need sum(nchar(funcderiv)) as width, and even then need to do some sort of
# paste(deparse(...),collapse='') to get a single string
    nrfunc <- paste(text='(',func,')/(',deparse(funcderiv, width=500),')', collapse='')
# first arg to outer()  will give rows
# Stupid Bug: I need to REVERSE zim to get proper axis orientation
    zstart<- outer(rev(zim*1i), zreal, "+")
    zindex <- 1:(length(zreal)*length(zim))
    zvec <- data.frame(zdata=as.vector(zstart), zindex=zindex,     itermap=rep(0,length(zindex)), badroot=rep(0,length(zindex)), rooterr=rep(0,length(zindex)) )

#initialize data.frame for zout.  
    zout=data.frame(zdata=rep(NA,length(zstart)), zindex=rep(NA,length(zindex)),     itermap=rep(0,length(zindex)), badroot=rep(0,length(zindex)), rooterr=rep(0,length(zindex)))
    # a value for rounding purposes later on; yes it works for  rootprec >1 
    logprec <-  -floor(log10(rootprec))
    newtparam <- function(zvar) {}
    body(newtparam)[2]  <- parse(text=paste('newz<-', nrfunc, collapse=''))
    body(newtparam)[3] <- parse(text=paste('return(invisible(newz))'))
    iter <- 1
    zold <- zvec  # save zvec so I can return original values
    zoutind <- 1 #initialize location to write solved values
    while (iter <= maxiter & length(zold$zdata)>0 ) {
        zold$rooterr <- newtparam(zold$zdata)
        zold$zdata <- zold$zdata - zold$rooterr
        rooterr <- abs(zold$rooterr)
        zold$badroot[!is.finite(rooterr)] <- 1
        zold$zdata[!is.finite(rooterr)] <- NA
# what if solvind = FFFFFFF? -- can't write 'nothing' to zout
        solvind <- (zold$badroot >0 | rooterr<rootprec)
            if( sum(solvind)>0 ) zout[zoutind:(zoutind-1+sum(solvind)),] <- zold[solvind,]
    #update zout index to next 'empty' row
        zoutind<-zoutind + sum(solvind)
# update the iter count for remaining elements:
        zold$itermap <- iter
# and reduce the size of the matrix being fed back to loop
        zold<-zold[!solvind,]
        iter <- iter +1
    # just wonder if a gc() call here would make any difference
# wow -- it sure does
        gc()
    }  # end of while
# Now, there may be some nonconverged values, so:
#  badroot[]  is set to 2  to distinguish from Inf/NaN locations
        if( zoutind < length(zindex) ) { # there are nonconverged values
#  fill the remaining rows, i.e. zout.index:length(zindex)
            zout[(zoutind:length(zindex)),] <- zold # all of it
            zold$badroot[] <- 2 # yes this is safe for length(badroot)==0
            zold$zdata[]<-NA #keeps nonconverged values from messing up results
            }
#  be sure to properly re-order everything...
    zout<-zout[order(zout$zindex),]
    zout$zdata <- complex(re=round(Re(zout$zdata),logprec), im=round(Im(zout$zdata),logprec) )
    rootvec <- factor(as.vector(zout$zdata), labels=c(1:length(unique(na.omit(as.vector(zout$zdata))))))
    #convert from character, too!
    rootIDmap<-matrix(as.numeric(rootvec), nr=length(zim))
# to colorize very simply:  
    if(drawplot) {
             colorvec<-rainbow(length(unique(as.vector(rootIDmap))))
        imagemat<-rootIDmap
        imagemat[,]<-colorvec[imagemat]  #now has color strings
        dev.new()
# all '...' arguments used to set up plot
        plot(range((zreal)),range((zim)), t='n',xlab='real',ylab='imaginary',... ) 
        rasterImage(imagemat, range(zreal)[1], range(zim)[1], range(zreal)[2], range(zim)[2], interp=F)     
        }

    outs <- list(rootIDmap=rootIDmap, zvec=zvec, zout=zout, nrfunc=nrfunc)
    return(invisible(outs))
}

Upvotes: 2

Josh O'Brien

Reputation: 162451

This first problem arises because readline reads in a text string, whereas what you need is an expression. You can use parse() to convert the text string to an expression:
```
f <-readline(prompt="Function? ")
sin(x)
f
# [1] "sin(x)"

f <- parse(text = f)
f
# expression(sin(x))

g <- D(f, "x")
g
# cos(x)
```
To pass in values for the arguments in the function call in the expression (whew!), you can eval() it in an environment containing the supplied values. Nicely, R will allow you to supply those values in a list supplied to the envir= argument of eval():
```
> eval(f, envir=list(x=0))
# [1] 0
```

Upvotes: 11

Henry

Reputation: 6784

Josh has answered your question

For part 2 you could have used

g <- expression( sin(x) )

g[[1]]
# sin(x)

f <- function(x){ eval( g[[1]] ) }

f(0)
# [1] 0
f(pi/6)
# [1] 0.5

Upvotes: 3

R type conversion expression() function()

Answers (3)

Related Questions