R - implicit external function command

Question

I would like to build an implicit function of x and y: 1 - 0.125 * y ^ 2 - x ^ 2 = 0.005 I know that it can be plotted as a contour plot, but it has problems with the "external" command in the following:

x<-seq(0.4,1.01,length=1000)
y<-seq(0,3,length=1000)
z<-outer(x,y,FUN="1-0.125*y^2-x^2=0.005")
contour(x,y,z,levels=0,drawpoints=FALSE)

I read the FAQ (7.17) regarding the "external" command and the need to vectorize the function, but I'm still in quandry.

Answer 1

I think you are a little confused about the meaning of "function". All operations (+, -, ^) in your function are vectorized so that everything works fine.

 x <- seq(0.4,1.01,length=1000) y <- seq(0,3,length=1000) z <- outer(x,y,function(x,y) 1-0.125*y^2-x^2-0.005) contour(x,y,z,levels=0,drawlabels=FALSE) 

Or if you need a small shortcut:

 library(emdbook) curve3d(1-0.125*y^2-x^2-0.005, xlim=c(0.4,1.01), ylim=c(0,3), n=c(100,100), sys3d="contour",drawlabels=FALSE,levels=0) 

This is actually slower because it uses the inner for loop, not outer() , so I set it to 100x100, not 1000x1000 (which is too large for this example anyway), but it will work on more complex examples cannot be easily vectorized ...

Answer 2

building an implicit equation as a contour graph is redundant. You essentially throw away 99.99% of the calculations you did.

It is best to find the y value for a given x that will make equation 0. Here is the code using uniroot in the R base.

R code using uniroot from the R base

 x = seq(0, 0.995, length = 100) # no real root above x = 0.995 root <- function(a) { uniroot(function(x,y) 1 - 0.125 * y^2 - x^2 - 0.005, c(0, 3), x = a )$root #only care about the root } y <- sapply(x, root) plot(x,y, type = "l") 

Ok c (0, 3) in the uniroot argument is the range of y values ​​where the root lies. therefore, for each given x value, uniroot will look for a y value from 0 to 3 for the root.

R code using fsolve from the pracma package

 library("pracma") x <- seq(0,0.995, length=100) fun <- function(y) c(1 - 0.125 * y^2 - x^2 - 0.005) y_sol <- fsolve(fun, rep(1, length(x)))$x plot(x,y_sol, type="l") 

fsolve takes a function whose root is the desired one and guesses the y values ​​for each given x value. Here we say that the values ​​of y lie near 1. We give it the value of assumption 1

uniroot wants the related range to look for root, fsolve needs to guess where the root may be.

These are faster ways to build implicit equations. Then you can use any graphics package, for example ggplot2 / rbokeh, to build a graph.

I have not done any tests, so I can’t say which one is faster. Although it does not matter for such a simple implicit function.