StopReadingThisUsername StopReadingThisUsername - 1 month ago 7
R Question

Integration in R giving "non-finite function value"

integral <- function(x) {2.393794315*((3320)*(x/2.24581)^5+(-1613880/11)*(x/2.24581)^4+(171163181/66)*(x/2.24581)^3+(-7563546913/330)*(x/2.24581)^2+(835541173981/8250)*(x/2.24581)+(-2953570085669/16500))*(((483793.161846485)*x^8+(-76823340.9717028)*x^7+(5337025908.822)*x^6+(-211866341077.587)*x^5+(5256530719898.47)*x^4+(-83466263852549.1)*x^3+(828318375700455)*x^2+(-4697211251008830)*x+(11653475160809900))^0.5)}

integrate(integral, lower=19.538547, upper=20.3245805)


This gave me


Error in integrate(integral, lower = 19.538547, upper = 20.3245805) :
non-finite function value


I am not sure what to do. It'd be great if somebody could input the integral into a software like Maple to see if it works - or, giving an answer to how to get around this error would be great, too :P

Thanks in advance!

EDIT:

The function that I am trying to revolve around the x-axis is:

y=2.393794315*((3320)(x/2.24581)^5+(-1613880/11)(x/2.24581)^4+(171163181/66)(x/2.24581)^3+(-7563546913/330)(x/2.24581)^2+(835541173981/8250)*(x/2.24581)+(-2953570085669/16500))

When plotted with limits {19.538547 < x < 20.3245805}:

enter image description here

As in the link provided by @StephenWade (in the comments), to find the surface area, I needed (f'(x))^2 + 1. Using Wolfram, Alpha, I computed this to be:

(f'(x))^2 + 1 = (483793.161846485)*x^8+(-76823340.9717028)*x^7+(5337025908.822)*x^6+(-211866341077.587)*x^5+(5256530719898.47)*x^4+(-83466263852549.1)*x^3+(828318375700455)*x^2+(-4697211251008830)*x+(11653475160809900)

Plotting this, without limits, gives:

enter image description here

However, computing this on R, I got an error:

enter image description here

My input into R was the following:

2.393794315*((3320)(x/2.24581)^5+(-1613880/11)(x/2.24581)^4+(171163181/66)(x/2.24581)^3+(-7563546913/330)(x/2.24581)^2+(835541173981/8250)(x/2.24581)+(-2953570085669/16500))(((483793.161846485)*x^8+(-76823340.9717028)*x^7+(5337025908.822)*x^6+(-211866341077.587)*x^5+(5256530719898.47)*x^4+(-83466263852549.1)*x^3+(828318375700455)*x^2+(-4697211251008830)*x+(11653475160809900))^0.5)

I just need a reasonable estimate of the surface area...Any ideas?

Answer

I would approach it using the polynom package, which can calculate the derivative for us, and then apply the formula from Wolfram.

revol_coef <- 2.393794315 * c(-2953570085669 / 16500,
                              835541173981 / (8250 * 2.24581),
                              -7563546913/ (330 * 2.24581^2),
                              171163181 / (66 * 2.24581^3),
                              -1613880 / (11 * 2.24581^4),
                              3320 / 2.24581^5)
y <- polynomial(revol_coef)
y_d <- deriv(y)

f <- function(x) {
    2 * pi * predict(y, x) * sqrt(1 + predict(y_d, x)^2)
}


integrate(f, lower = 19.538547, upper = 20.3245805)

The output I get is

-45.71118 with absolute error < 0.0017

The answer is negative, due to the surface being specified in the lower-right quadrant of the x-y plane (as per your plot). To fix this, take the negative of the number supplied to get a surface area.

Translating formulas from Wolfram Alpha, manipulating them, and then putting them into R can easily involve making a mistake.

In this case, I would recommend doing as much of the work within R, as the tools/functions to do the necessary calculations and simplifications are readily available.

Comments