I have been using R for about two years, with no previous coding background. So, I feel like the title says, “always learn and never know”. This time, I decided to use R to study a simple, non-statistical problem that came up some time ago.

Suppose the exponential function 2^x and the parabola x^2. One can readily show, e.g with the Mean Value Theorem, that the two curves intersect in exactly three points, at x=2, at x=4 and at a third position in the interval (-1,0). We are interested in this third negative root.

png(filename = "plot1.png", width = 400, height = 300)
x <- seq(-2,5,0.1)
mt = "Intersection plot"
plot(x, 2^x, type="l", xlim=range(x), main=mt, ylab="y", ylim=c(0,30), lwd=2)
lines(x, x^2, type="l", cex=1.5, lty=2, lwd=2)
exp1 <- expression(paste(italic(y)*"="*italic(x)^2))
exp2 <- expression(paste(italic(y)*"="*2^italic(x)))
legend(-2, 30, c(exp1, exp2), cex=1.5, col=rep("black",2), lty=c(1,2), lwd=rep(2,2))
points(c(-0.7666647, 2, 4), c(0.5877748, 4, 16), col="red", pch=19, cex=2)
dev.off()

Now the question is: given an exponential base a>=1, there is always a negative root root(a). How can we plot this function root?

When it comes to custom functions, R offers remarkable capabilities. Here, we can use a standard method, such as Newton-Raphson. After some plain calculus, we reach to an expression that provides the next approximation of the root:

(a^x*(x*log(a)-1)-x^2)/(a^x*log(a)-2*x)

Here’s how to build this in R. First, we define a new function, call it root, and we also define the variable name a. After the variable a has been set, we can construct the approximation formula, given above. Then, we create an empty vector to store the sequence of root approximations and we give the first value. Note that this value has to be a “good” one. A number of 10 steps can yield a quite satisfactory approximation. The last step is to indicate the value for the defined root function, with the return command. So, we have a function defined into a function.

root <- function(a) {
newtonsq <- function(x) {(a^x*(x*log(a)-1)-x^2)/(a^x*log(a)-2*x)}
rootsq <- c()
rootsq[1] <- -0.5
for (j in 1:9) { rootsq[j+1] <- newtonsq(rootsq[j]) }
return(rootsq[10])
}

Now try our new function. You can notice that it is strictly increasing, assuming values in [-1,0). However, the increasing rate is very slow, so we need a scaled x-axis. The result has an hyperbolic shape.

png(filename = "plot2.png", width = 400, height = 300)
a <- c(1:10, 10^(1:100))
b <- sapply(a, root)
plot(log(a), b, type="l", xlim=c(0,100), ylim=c(-1,0), main="Negative root", ylab="root", cex=1.5, lwd=2)
dev.off()

The equation is related to the Lambert function. An interesting discussion about the equation can be found here. Finally, there is a built-in function in R, called uniroot, that can be used in place of the custom function defined herein.