R type conversion выражение() function()

Question

Я пытался написать программу в R, которая реализует метод Ньютона. Я был в основном успешным, но есть две маленькие коряги, которые меня беспокоят. Вот мой код:

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) 
} 

1. Я не могу получить R распознавать функции f и f. если я пытаюсь использовать readline() для запроса ввода пользователя. Я получаю сообщение об ошибке «Ошибка в Newton(): не удалось найти функцию« f ».« Однако, если я прокомментирую строки чтения (как указано выше), заранее определите f и f., тогда все будет хорошо.

2. Я пытался заставить R вычислить производную от функции. Проблема в том, что объект класса, с которым R может принимать символические производные, равен expression(), но я хочу взять производную от function() и дать мне function(). Короче говоря, у меня возникают проблемы с преобразованием типов между expression() и function().

У меня есть уродливые, но эффективное решение для перехода от function() к expression(). При заданной функции f D(body(f)[[2]],"x") даст производную от f. Однако этот вывод является expression(), и я не смог вернуть его обратно в function(). Нужно ли использовать eval() или что-то в этом роде? Я пробовал подмножество, но безрезультатно. Например:

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

когда то, что я хочу, е (0) = 0, так как грех (0) = 0.

EDIT: Спасибо всем! Вот мой новый код:

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) 
} 

Answer 1

1. Это первая проблема возникает потому, что readline читает в текстовой строке, в то время как то, что вам нужно, это выражение. Вы можете использовать parse() преобразовать текстовую строку в выражение:

   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) 
   

2. Чтобы пройти в значениях для аргументов при вызове функции в выражении (! Гмм), вы можете eval() его в среде, содержащей входящий в комплект поставки значения. Ну, R позволит поставить эти значения в списке, подводимый к envir= аргументу eval():

   > eval(f, envir=list(x=0)) 
   # [1] 0 
   

Answer 2

Джош ответил на ваш вопрос

В части 2 вы могли бы использовать

g <- expression(sin(x)) 

g[[1]] 
# sin(x) 

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

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

Answer 3

BTW, недавно написав игрушку, которая вычисляет фрактальные узоры, основанные на корневой конвергенции метода Ньютона в комплексной плоскости, я могу порекомендовать вам бросить так me code like the following (где список аргументов основной функции включает в себя «func» и «varname»).

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

Если вы более осторожны, вы могли бы включать в качестве аргумента «funcderiv», и оберните мой код в

if(missing(funcderiv)){blah blah} 

Ahh, почему нет: вот моя полная функция для всех, чтобы использовать и пользуются :-)

# 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)) 
}