Peng Peng
Peng Peng

Reputation: 1365

Solve equation, strange result

here is my code to solve equation:

fx=function(x){ x^3-x-3}  
solve=function(a,b,eps){
    if(abs(fx(a))<0.00001)  return(list(root=a,fun=fx(a)))    
    else if(abs(fx(b))<0.00001)  return(list(root=b,fun=fx(b)))  
    else if (fx(a)*fx(b)>0)  return(list(root="failed to find"))  
    if (a>b){
        c<-a
        a<-b
        a<-b}       
    while( b-a>eps ){ 
        x=(a+b)/2
        if (fx(x)==0) {return(list(root=x,fun=fx(x))) } 
        else if (fx(a)*fx(x)<0) {b=x }           
        else  {a=x}}
    myroot=(a+b)/2
    return(list(root=myroot,value=fx(myroot)))
}

> solve(1,3,1e-8)
$root
[1] 1.6717

$value
[1] 2.674228e-08

> fx(1.6717)
[1] 8.73813e-07

Why fx(1.6717) != $value, I want to know the reason
8.73813e-07!=2.674228e-08?

how can i revise:return(list(root=myroot,value=fx(myroot)))
to make my root more digits ?

Upvotes: 1

Views: 247

Answers (2)

Andrie
Andrie

Reputation: 179398

When R prints a value, it uses by default digits=3, i.e. printing 3 significant digits. This means you made an error of interpretation when looking at your results.

Try this:

 x <- solve(1,3,1e-8)

print(x[[1]], digits=9)
[1] 1.67169989

Now substitute the actual returned value into your function:

fx(x[[1]])
[1] 2.674228e-08

Now the values match.

In summary, you have made a rounding error when interpreting the printed results of your function.


You can trace this behaviour in the R help files as follows:

?print

will point you to

?print.default

Which has this to say about the digits argument:

digits: a non-null value for digits specifies the minimum number of significant digits to be printed in values. The default, NULL, uses getOption(digits). (For the interpretation for complex numbers see signif.) Non-integer values will be rounded down, and only values greater than or equal to 1 and no greater than 22 are accepted.

Upvotes: 2

Alan
Alan

Reputation: 3203

Try this and look at the print() of a and b.

fx=function(x){ x^3-x-3}  
solve=function(a,b,eps){
    if(abs(fx(a))<0.00001)  return(list(root=a,fun=fx(a)))    
    else if(abs(fx(b))<0.00001)  return(list(root=b,fun=fx(b)))  
    else if (fx(a)*fx(b)>0)  return(list(root="failed to find"))  
    if (a>b){
        c<-a
        a<-b
        a<-b}       
    while( b-a>eps ){ 
        x=(a+b)/2
        if (fx(x)==0) {return(list(root=x,fun=fx(x))) } 
        else if (fx(a)*fx(x)<0) {b=x }           
        else  {a=x}}
    myroot=(a+b)/2
    print(a,digits=20)
    print(b,digits=20)
    return(list(root=myroot,value=fx(myroot)))
}

There is a round.

Upvotes: 0

Related Questions