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:

Capture d’écran_2023-01-15_12-19-00

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:

Capture d’écran_2023-01-15_12-19-37

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

This is not a Julia related question. Better ask here: https://groups.google.com/g/sympy?pli=1

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
       Base.adjoint(f::Function) = Base.Fix1(ForwardDiff.derivative, f)

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

julia> 2π*quadgk(g, 1, 3)[1]
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)
dq = 2π*quadgk(h, 1, 3)[1]

println("(Area, error) = ", d)

println("Area with QuadGK = ", dq)

(Area, error) = (325.4437377507612, 5.56797585815616e-9)
Area with QuadGK = 325.44373775078094