Plotting implicit functions in R

So in prepping for my latest manuscript on population dynamics I have been creating all the necessary figures.  One of them I considered was a 2-d surface plot of a modified Ricker equation showing the transitions from extinction stability, and stability to limit cycles.  Inconveniently though the only way to do this is with an implicit function.  Since becoming a fully ggplot2 convert, I wanted to be able to do it in ggplot2. Let’s start with a simple definition so we’re all clear. Here’s an explicit function. \[y=2^2+x^2\] It’s explicit because y is stated explicitly in terms of x, i.e. \[f(x)=2^2+x^2\] Yet sometimes its not possible to define y explicitly in terms of x. Implicit functions are formally defined as equations that satisfy the condition: \[R(x,y)=0\] A good example is the unit circle. The unit circle can be found in the equation: \[f(x,y)=x^2+y^2\] And unit circle is created when: \[f(x,y)=1\] Therefore we can draw the unit circle using the equation: \[x^2+y^2=1\] \[x^2+y^2-1=0\] This is the implicit function for the unit circle. You can solve it explicitly for y, but it requires two functions. \[y=\sqrt{1-x^2}, y=-\sqrt{1-x^2}\] So, if you have an implicit function you want to plot, how could you do it in R? The easy answer is to numerically solve your equation over a specific range, but then how can you plot it?

my.fun <- function(x,y){y^2+x^2 -1}
x<-seq(-1.5,1.5,length=1000)
y<-seq(-1.5,1.5,length=1000)
z<-outer(x,y,my.fun)
contour(x,y,z,level=0)

That will produce a graph that looks like the figure on the left, whereas the figure on the right was created with the function I wrote called imp.solve().

my.fun <- function(x,y){y^2+x^2 -1}
to.plot <- imp.solve(x,y,my.fun)
#Plot with ggplot
ggplot(to.plot,aes(x=x,y=y))+geom_point()+xlim(range(x))+ylim(range(y))

The first example produces a quick and dirty plot that is all fine and good, but its a surface plot. That means you’re stuck with the “0” on the line denoting the surface is only drawing areas where 0 is. On top of that you lose lots of other control over what you can draw and what if you you want to use ggplot2, or some other package, or you actually want to know what the numerical answer is? My function parses
 the output from the outer() function and returns a dataframe you can easily plot in ggplot2. It can handle most polynomial functions, although I’ll leave you in charge of ordering your own output for line plots in scenarios with multiple x’s at the same y point. Here’s another example with a cubic function. \[x^3+y^3-1=0\]







So why does this matter in ecology? Often you might want to create a surface plot for transitions between stable states and unstable states in with a discrete time equation. Here’s an example from my own work. I’m trying to recreate a figure such as this from a paper by Aviles 1999. Here is the figure at right. The function is a modified Ricker equation: \[N_{t+1}={N_t}^{(1+\gamma)}e^{(-c+N_t)}e^r\] Here gamma is the allee effect term, c is the strength of density dependence and r is the population growth rate. Both gamma and r are bifurcation points and we can use an implicit function to draw the boundaries between say stability and extinction at a constant c. So while implicit functions aren’t the most common thing to come across in ecology, they do arise with some mathematical models. In this example the transition between stability and extinction is given by the equation: \[ce^{(\gamma – r)}-\gamma=0\] Holding c constant at .01, the following code will generate our figure.

##########First set up the parameters to vary over########
y<-seq(0.01,.9,length=1000)
x<-seq(-2,8,length=1000)
###################Define the function############
my.fun <- function(x,y){.01*exp((y-x)/y)-y }
to.plot <- imp.solve(x,y,my.fun)

ggplot(to.plot,aes(x=x,y=y))+geom_line()

This is the same as the the first line in the above figure from the paper. You can find the full code for this function over at my github page as Imp.solve.R. So while a previously quick and dirty solution existed in R, here is an improved version with greater plotting flexibility. Happy implicit function graphing!