Fast and Sloppy Root Finding

Disclaimer: While the approach to root finding mentioned below is both slow and imprecise, it's also a cheap and incredibly handy approach when all you need to do is get "close enough". If you like R quick and dirty (hey now, get your mind out of the gutter...) this is totally the root finding method for you!

I just read a post on Root Finding (original here) by way of R-bloggers.com which was timely given that only yesterday I'd needed to do some root finding in R to make a figure for a manuscript -- something like the following image.

The blog post prompted me to mention here how I did my root finding for two reasons:

1. Precision and computation time sometimes don't matter all that much; and
2. The way I did my root finding was way easier to implement than convergence-based methods described in the post above.

So here's what I was aiming for, and how I implemented it in R.

The Task: Suppose you're plotting two curves (say, y=f(x) and y=g(x)) and would like to indicate their intersection with an empty circle (i.e. pch=21). In my case, the intersection of these two curves was equilibrium point for a dynamic model, and I wanted to indicate it as such.

If you can find their intersection mathematically (i.e. set f(x)=g(x) and solve for x) then awesome -- do that if you can.  But if for some reason you can't, and you know a single root exists in some interval a≤x≤b, you can find that root quickly using some straightforward vector tricks.

The Solution: Lets use the example of finding the intersection of f(x) = x/(1+x) and g(x) = (5-x)/5 over the interval (0,5).

Step 1:  Define an x vector full of values between a and b.  The smaller the step size (or, the longer the list) the better.

> x = seq(0, 5, length=5000);

Step 2: Compute the square of the difference of your two functions over that interval using x.  This is as simple as the line of code...

> fgdiff = ( x/(1+x) - (5-x)/5 )^2;

Step 3: Using the which() function, we can pick out the index for the smallest value in our list of squared differences... Once we know this index (call it j) we know the intersection occurs at, or very near, the x value x[j], and we're basically done! 

> j = which(fgdiff==min(fgdiff))[1]; 
> j; x[j];  ## show the value of j, x[j]

For a closer look at what's going on with that which() statement, check out the help for which() and following example.

## A Closer look at the which(fgdiff==min(fgdiff))
> ?which
> which(c(F,F,T,T,F))
[1] 3 4
> which(c(F,F,T,T,F))[1]
[1] 3
> xample = c(5:1,2:5); xample
[1] 5 4 3 2 1 2 3 4 5
> min(xample)
[1] 1
> xample==min(xample)
[1] FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE
> which(xample==min(xample))
[1] 5
> which(xample==min(xample))[1]
[1] 5

Step 4: Since both functions are (approximately) equal at this x value, it only remains to decide whether you want to indicate the point of intersection using (x, f(x)) or (x, g(x)).

> points(x[j], x[j]/(1+x[j]), pch=19, cex=2)

All done!

If you'd like to tinker with this example, here's the code to produce the image above.

win.graph(5,4)   ## Open a new window, 5in x 4in.

## Plot the two functions
curve(x/(1+x), from=0, to=5, lwd=3, ylab="f(x), g(x)", 
                             ylim=c(0,1), cex.lab=1.5)
curve((5-x)/5, from=0, to=5, lwd=3, lty=2, add=T)

## Find their intersection the quick n' dirty way
x=seq(0,5,length=1000)
fgdiff = (x/(1+x) - (5-x)/5)^2;
j = which(fgdiff==min(fgdiff))[1]; # Only 1st, in case of 2 mins
j; x[j]

## Plot an empty circle at the intersection, and label curves.
points(x[j], x[j]/(1+x[j]), pch=21, col="black", bg="white", lwd=3, cex=1.7)
text(4.5, 0.65,  expression(f(x)==frac(x,1+x)), cex=1)
text(3,   0.15, expression(g(x)==frac(5-x,5)), cex=1)