Always learn and never know

[This article was first published on infoRύχος » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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.

To leave a comment for the author, please follow the link and comment on their blog: infoRύχος » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)