Area of Surface of Revolution Integral too Hard to be Computed By JULIA' SymPy and Python' SymPy

Hi all,

I want to get the numerical answer of the area integral.

the curve is y =( x^6 + 2)/(8x^2) with 1 \le x \le 3, revolved about x-axis.

the formula is:
A = 2 \pi * \int_{a}^{b} f(x) * \sqrt{1 + f'(x)} dx

Julia code:

using SymPy

x = symbols("x")

f(x) = (x^6 + 2)/(8x^2)
g(x) = sqrt(1 + diff(f(x),x))

h = 2pi*integrate(((x^6 + 2)/(8x^2))*sqrt(1 + diff(f(x),x)), (x, 1, 3))
d = simplify(h)


the result is in this form:

with Python 3.9:

import sympy as sy

x = sy.Symbol("x")

def f(x):
return ((x**6) + 2)/ (8*x ** 2)

def fd(x):
return sy.simplify(sy.diff(f(x), x))

def f2(x):
return sy.sqrt((1 + (fd(x)**2)))

def vx(x):
return 2*np.pi*(f(x)*((1 + (fd(x) ** 2))**(1/2)))

vx = sy.simplify(sy.integrate(vx(x), (x, 1, 3)))



the result is like this:

is this integral of area of surface too hard to be computed?

1 Like

That integral looks hard and I don’t see an obvious analytical solution. If you want a numerical answer, then numerically evaluate it using QuadGK or something similar.

1 Like

Thanks for the mailing list

Let’s try numerical integration.

It takes three lines in Julia to make a trapezoidal integration function, out of which one line is end:

integrate(f, range) = let r=range, a=first(r), b=last(r), dx=step(r)
(0.5f(a) + sum(f, r[2:end-1]) + 0.5f(b))dx
end


to check that it works:

julia> integrate(abs2, range(0,1,10^8))
0.3333333333333333


For the calculation of f^\prime(x) I’m going to use ForwardDiff, which uses dual numbers to automatically differentiate the function.

julia> using ForwardDiff

julia> f(x) = (x^6+2)/8x^2
f (generic function with 1 method)

julia> g(x) = f(x)*√(1+f'(x))
g (generic function with 1 method)


Ok let’s try integrating your function.

julia> 2π*integrate(g, range(1, 3, 10^8))
116.28129729349027


Let’s confirm against QuadGK, which uses magic to run much faster for the same precision (and even offers information about its numerical precision!)

julia> using QuadGK

116.2812972934902


Julia saves the day yet again (when does it not?).

2 Likes

Wow amazing,

I agree Julia saves the day. Only a matter of creating plot of solid of revolution like Python that is still I haven’t figure out with Julia Plots.

Julia is fast in calculating this integral.

The result with using QuadGK with including error tolerance is different:

using SymPy, QuadGK

x = symbols("x")

f(x) = (x^6 + 2)/(8x^2)
# simplify(diff(f(x),x))
# we cannot use diff inside the fd(x) since QuadGK can't comprehend diff
fd(x) = (x^6 - 2)/(2x^3)
g(x) = sqrt(1 + (fd(x))^2)
#g(x) = f(x)*√(1+(fd(x))^2)

h(x) = f(x)*g(x)
Area(x) = 2pi*h(x)

d = quadgk(x -> Area(x), 1, 3, rtol=1e-10)