[R] how to see what's wrong with a self written function?

[Petr Savicky]

On Tue, Dec 21, 2010 at 11:39:31AM -0800, casperyc wrote:
> 
> Hi all,
> 
> I am writing a simple function to implement regularfalsi (secant) method.
> 
> ###################################################
> regulafalsi=function(f,x0,x1){
> 	x=c()
> 	x[1]=x1
> 	i=1
> 	while ( f(x[i])!=0 ) {
> 		i=i+1
> 		if (i==2) {
> 			x[2]=x[1]-f(x[1])*(x[1]-x0)/(f(x[1])-f(x0))
> 		} else {
> 			x[i]=x[i-1]-f(x[i-1])*(x[i-1]-x[i-2])/(f(x[i-1])-f(x[i-2]))
> 		}
> 	}
> 	x[i]
> }
> ###################################################
> 
> These work fine,
> regulafalsi(function(x) x^(1/2)+3*log(x)-5,1,10)
> regulafalsi(function(x) x^(1/2)+3*log(x)-5,10,1)
> 
> For all x>0, the function is strictly increasing.
> 
> Then
> 
> regulafalsi(function(x) x^(1/2)+3*log(x)-5,1,100)
> 
> Error in while (f(x[i]) != 0) { : missing value where TRUE/FALSE needed
> In addition: Warning message:
> In log(x) : NaNs produced
> 
> I dont know what happened there, is there a way to find the value for
> f(x[i])
> that R can't determine TRUE/FALSE?

The previous replies imply that at some point, your code evaluates the
function outside the required interval. The reason for this may be seen
using graphics. If you add plot commands to your code, for example

  regulafalsi=function(f,x0,x1){
          z <- seq(x0, x1, lenth=1000)
          plot(z, f(z), type="l")
          abline(h=0)
          x=c()
          x[1]=x1
          i=1
          while ( f(x[i])!=0 ) {
                  i=i+1
                  if (i==2) {
                          x[2]=x[1]-f(x[1])*(x[1]-x0)/(f(x[1])-f(x0))
                          lines(c(x0, x[i-1]), c(f(x0), f(x[i-1])), col=2)
                  } else {
                          x[i]=x[i-1]-f(x[i-1])*(x[i-1]-x[i-2])/(f(x[i-1])-f(x[i-2]))
                          lines(c(x[i-2], x[i-1]), c(f(x[i-2]), f(x[i-1])), col=2)
                  }
                  points(x[i], 0)
                  aux <- readline("press Enter to continue")
          }
          x[i]
  }
  regulafalsi(function(x) x^(1/2)+3*log(x)-5,1,10)
  regulafalsi(function(x) x^(1/2)+3*log(x)-5,1,100)

then it may be seen that for [1, 10], we get convergence, while for [1, 100]
the secant line in the second iteration does not intersect x-axis inside the
required interval.

Petr Savicky.