Solve equations in R, similar to the Excel solver options function

I have a question about the possibility of solving functions in R, and am doing the same using excel.

However, I want to do it with R to show that R is better for my colleagues :)

Here's the equation:

f0<-1e-9
t_pw<-30e-9
a<-30.7397582453682
c<-6.60935546184612

P<-1-exp((-t_pw)*f0*exp(-a*(1-b/c)^2))

      

I want to find a value b

for P<-0.5

. In Excel, we can do this by selecting the P values ​​column and setting it to 0.5 and then using the solver options function.

I don't know which method is better? Or any other way to do it?

Thankx.

+3


source to share


2 answers


I have a strong suspicion that your equation should include -t_pw/f0

, not -t_pw*f0

, and t_pw

should be 3.0e-9

, not 30e-9

.

 Pfun <- function(b,f0=1e-9,t_pw=3.0e-9,
                  a=30.7397582453682,
                  c=6.60935546184612) {
               1-exp((-t_pw)/f0*exp(-a*(1-b/c)^2))
           }

      

Then @Lyzander's suggestion works uniroot()

great:

 u1 <- uniroot(function(x) Pfun(x)-0.5,c(6,10))

      



The estimate here is 8.05.

 par(las=1,bty="l")
 curve(Pfun,from=0,to=10,xname="b")
 abline(h=0.5,lty=2)
 abline(v=u1$root,lty=3)

      

enter image description here

+2


source


If you want to solve an equation, the simplest is to use uniroot

which is in base-R.

f0<-1e-9
t_pw<-30e-9
a<-30.7397582453682
c<-6.60935546184612

func <- function(b) {
    1-exp((-t_pw)*f0*exp(-a*(1-b/c)^2)) - 0.5
}

#interval is the range of values of b to look for a solution
#it can be -Inf, Inf
> uniroot(func, interval=c(-1000, 1000), extendInt='yes')
Error in uniroot(func, interval = c(-1000, 1000), extendInt = "yes") : 
  no sign change found in 1000 iterations

      

As you can see above, my function unitroot

fails. This is because there is no single solution to your equation that is easy to see. exp(-0.0000000000030 * <positive number between 0-1>)

practically (very close) 1, so your equation becomes 1 - 1 - 0.5 = 0

that doesn't hold. You can see the same with the plot:

curve(func) #same result for curve(func, from=-1000, to=1000)

      

enter image description here



In this function, the result will be -0.5 for any b.

So, one way to do it quickly uniroot

, but probably for a different equation.

And a working example:

myfunc2 <- function(x) x - 2 

> uniroot(myfunc2, interval=c(0,10))
$root
[1] 2

$f.root
[1] 0

$iter
[1] 1

$init.it
[1] NA

$estim.prec
[1] 8

      

+3


source







All Articles