Skip to content

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

5 messages · casperyc, Duncan Murdoch, jim holtman +1 more

casperyc
# 
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?

Thanks!

casper
Duncan Murdoch
#
On 21/12/2010 2:39 PM, casperyc wrote:
The easiest is to just use regular old-fashioned debugging methods, i.e. 
insert print() or cat() statements into your function.  You could also 
try debug(regulafalsi) and single step through it to see where things go 
wrong.  (An obvious guess is that one of the values being passed to f is 
negative, but you'll have to figure out why that happened and what to do 
about it.)

Duncan Murdoch
jim holtman
# 
Here is what I get when I have:

options(error=utils::recover)

I always run with the option so that on an error, I get dumped in the
browser to see what is happening.  It appears that 'i == 3' when the
error occurs and you can also see the values of 'x':
+        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]
+ }
[1] 2.978429
Error in while (f(x[i]) != 0) { : missing value where TRUE/FALSE needed
In addition: Warning message:
In log(x) : NaNs produced

Enter a frame number, or 0 to exit

1: regulafalsi(function(x) x^(1/2) + 3 * log(x) - 5, 1, 100)

Selection: 1
Called from: top level
Browse[1]> i
[1] 3
Browse[1]> x
[1] 100.00000  18.35661 -42.22301
Browse[1]>
On Tue, Dec 21, 2010 at 2:39 PM, casperyc <casperyc at hotmail.co.uk> wrote:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Petr Savicky
    

    

      
      #
    

  


  

      
On Tue, Dec 21, 2010 at 11:39:31AM -0800, casperyc wrote:

            

      
      
      
    

        
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.

  
  


  


      

    

    
8 days later

    

      
      

        


  

        

      casperyc
    

    

      
      #
    

  


  

      
Hi Petr Savicky,

It is very interesting and a good way to plot,
so that helps me to visualized the method,
also easier to see what's wrong.

The "readline" is new to me,
and I should learn how to use it in the future.

Thanks!

casper