[R] How to solve difficult equations?

[(Ted Harding)]

On 25-Apr-07 07:15:55, francogrex wrote:
> 
> This below is not solvable with uniroot to find "a":
> fn=function(a){
> b=(0.7/a)-a
> (1/(a+b+1))-0.0025
> }
> uniroot(fn,c(-500,500))  gives
> "Error in uniroot(fn, c(-500, 500)) : f() values at end points
> not of opposite sign"
> 
> I read R-help posts and someone wrote a function:
> http://finzi.psych.upenn.edu/R/Rhelp02a/archive/92407.html
> but it is not very precise. Is there any '"standard" function
> in R that can solve this? thanks.

Two answers: Yes, and No.


First, "No":

Let alpha denote 0.0025, and beta 0.7 (in your function "fn").
Then

  fna <- function(alpha,beta){ beta*alpha/(1 - alpha) }

solves it. But this is not a standard R function.


Second, "Yes":

and the standard R function is uniroot(). But you can only apply
it usefully if you first study the behaviour of your function fn(),
in rather careful detail.

Over your range (-500,500):

  a<-10*(-50:50)
  plot(a,fn(a),pch="+")

Clearly something extreme happens just to the left of a=0. So:

  a <- 0.025*(-100:0)
  plot(a,fn(a),pch="+")

and so for this set of values of 'a' the previous behaviour
cannot be seen. So:

  a <- 0.01*(-100:100)+0.001
  plot(a,fn(a),pch="+")

so the function goes very negative somewhere around a = -0.7.
But

  fn(500)
  [1] 0.996102

so it is positive for a=500. Now find (inspired by the latest
plot):

  a[which(fn(a) < (-100))]
  [1] -0.699

and now you can use uniroot:

  uniroot(fn,c(-0.699,500))
  $root
  [1] 0.001771128
  $f.root
  [1] 2.379763e-05
  $iter
  [1] 16
  $estim.prec
  [1] 6.103516e-05

and, if that doesn't look precise enough:

  uniroot(fn,c(-0.699,500),tol=1e-10)
  $root
  [1] 0.001754386
  $f.root
  [1] 1.354602e-14
  $iter
  [1] 18
  $estim.prec
  [1] 5e-11


Now compare with the function fna() that solves it directly:

  fna(0.0025,0.7)
  [1] 0.001754386

(so in fact it was worth increasing the precision for uniroot).


But the lesson to be drawn from all this is that for functions
like fn(), which have singularities (here at a = -0.7), the
blind application of root-finding functions may not work, since
they are not set up to explore the function is the kind of way
illustrated above. While there are procedures in the numerical
analysis world to handle this kind of thing, they tend to be
written for particular classes of function, and again you will
have to do a bit of exploration to find out which function to use.

And (while someone more knowledgeable may well disagree with me)
I suspect that these are not "standard" R funnctions.

Hoping this is helpful,
Ted.


--------------------------------------------------------------------
E-Mail: (Ted Harding) <ted.harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 25-Apr-07                                       Time: 09:31:29
------------------------------ XFMail ------------------------------