2012-01-13 8 views
9

私は、ニュートンの方法を実装するRでプログラムを書くことを試みてきました。私はほとんど成功しましたが、私に迷惑をかけている2つの小さな障害があります。ここに私のコードは次のとおりです。R型変換式()function()

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. 私は、ユーザー入力を要求するようにreadline()を使用しようとするとRは機能ff.を認識することができません。しかし、私が上記のようにreadlineをコメントアウトして、ff.をあらかじめ定義しておけば、すべて正常に動作します。

  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) 

私がしたいことはfが(0)= 0以来の罪(0)= 0

編集:おかげですべての!ここに私の新しいコードです:

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

答えて

7
  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()することができます値。

    > eval(f, envir=list(x=0)) 
    # [1] 0 
    
あなたが使用している可能性がジョシュは、パート2の場合は、あなたの質問

に答えた

+0

:-)お楽しみいただき、ありがとうございます! readline()のより良い代替手段はありますか?私はasseを試みましたが、parse()についてはわかりませんでした。また、body()以外のfunction()からexpression()に行くためのより良いオプションがありますか? – Quasicoherent

+0

'readline()'はユーザからの型付き入力を行うための正しい関数です。また、 'parse()'よりも良い方法はありません。しかし、あなたが 'function()'から 'expression()'に行くことによってあなたが意味することを確信していません。 'function()'は関数呼び出しか関数定義ですか? –

+1

私は 'function()'クラスにあるものを意味し、特に数学的関数です。例えば、関数f -function(x){3x * sin(x/3)+ 2 * cos(4 * x)} 'を持っていて、その微分を求めたいとします。私がこれを行うために知っている唯一の方法は 'body()'を使って式を取得し、 'D(body(f)[[2]]、" x ")'を使うことです。 – Quasicoherent

1

:うまく、Rは、あなたが eval()envir=引数に供給されたリストにそれらの値を供給することができます
g <- expression(sin(x)) 

g[[1]] 
# sin(x) 

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

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

ありがとうございました!もう一つのしわ:表現変数(例ではg)を取り除きたいのですが? {([[1]] F)はeval} それが評価」であるため、円形の依存関係を作成する((x)関数 - 式(SIN(X)) F < - 実際に、私は F <を実行したいですf [[1]]) 'f [[1]]はfで実際に動作するもの(すなわち' sin(x) ')に置き換えられません。 – Quasicoherent

+0

私はもはやあなたが求めているものは明確ではない。おそらくあなたはgを持たないこのようなものを探しているでしょう: 'f < - function(x){eval(expression(sin(x))[[1]])}; f(pi/6) 'あるいは多分' x < - pi/4;関数ではないもの; eval(expression(sin(x))) ' – Henry

+0

まさに - gを消したいです。私はgの内容を抽出してからgを取り除きたいと思っています。私の例を使って: '> g <-expression(sin(x)) > f <-function(x){eval(g [[1]]) > f function(x){eval(g [ [1]]) '' 私が本当に望むのは '> f function(x){sin(x)}' つまり、fのg依存関係を取り除きたいと思います。それは理にかなっていますか? – Quasicoherent

2

最近、複素平面上のニュートン法の根収束に基づいてフラクタルパターンを計算するおもちゃを書いて、私のコードは次のようになります(main関数の引数リストには "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} 

ああに私のコードをラップすることができ、なぜない:すべてを使用するためにここに私の完全な機能だと

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

コードを送信していただきありがとうございます。私はまだそれを読んでいませんが、コードをエミュレートしてプログラムを改善できると確信しています。 – Quasicoherent