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

[Freya_the_Goddess]

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?

[ufechner7]

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

[Ininterrompue]

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.

[Freya_the_Goddess]

Thanks for the mailing list

[uniment]

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?).

[Freya_the_Goddess]

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.

[Freya_the_Goddess]

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

[oscarbenjamin]

@Freya_the_Goddess

(Area, error) = (325.4437377507612, 5.56797585815616e-9)

The area of revolution formula should have (1 + f'(x)^2)^{\frac{1}{2}}. The answer 116 comes from missing off the square inside the brackets whereas the correct answer is closer to 325 although that result seems to have a significant approximation error.

I’m not sure how exactly to do this from Julia but using SymPy from Python it would be:

from sympy import *

x = symbols('x')
f = (x**6 + 2) / (8*x ** 2)
A = 2*pi*Integral(f*sqrt(1 + f.diff(x)**2), (x, 1, 3))
print(A.evalf())

This gives 326.919561445782. This method of computing the integral numerically is slow but accurate. In particular it should be possible to trust that it is accurate to the number of digits that are shown in the output:

>>> A.evalf(50)
326.91956144578231119755087750201147914688815885596

It is also possible to compute this symbolically but SymPy’s integrate function needs the integrand to be simplified a little. First declare x to be positive and secondly apply a simplification inside the square root with sqf (square free factorisation):

from sympy import *

x = symbols('x', positive=True)
f = ((x**6) + 2)/ (8*x ** 2)
A = 2*pi*Integral(f*sqrt(sqf(1 + f.diff(x)**2)), (x, 1, 3))
print(A.doit())

Now the integral has simplified to

2 \pi \int\limits_{1}^{3} \frac{\left(x^{6} + 1\right) \left(x^{6} + 2\right)}{16 x^{5}}\, dx

and symbolic evaluation with doit gives 8429*pi/81 which agrees with the numerical result for the integral.

[Freya_the_Goddess]

Thanks for the reply, there is never too late for a reply. I tried this and works with PyCall in Julia.

this is the Julia code, but you have to have PyCall and the ENV set up to Python inside Conda installed in Julia… a bit of work but very rewarding to be able to test Python code in Julia:

# Calculate the surface area of y = (x**6 + 2) / (8*x ** 2)
# revolved about the x-axis
# This method of computing the integral numerically is slow but accurate. 
# In particular it should be possible to trust that it is accurate 
# to the number of digits that are shown in the output:
# A.evalf(50)

# First declare x to be positive and secondly apply a simplification 
# inside the square root with sqf (square free factorisation):
using PyCall

ENV["PYTHON"] = "/home/browni/.julia/conda/3/x86_64/bin/python3"
# the path is from the command 'which python3'

py"""
from sympy import *

x = symbols('x')
f = (x**6 + 2) / (8*x ** 2)
A = 2*pi*Integral(f*sqrt(1 + f.diff(x)**2), (x, 1, 3))
# A.evalf(50)

x = symbols('x', positive=True)
f = ((x**6) + 2)/ (8*x ** 2)
A = 2*pi*Integral(f*sqrt(sqf(1 + f.diff(x)**2)), (x, 1, 3))
print('Area =', A.doit())
print('=',A.evalf())

"""

[j_verzani]

This can be done with SymPy.jl, leaving the details of PyCall and the slightly different set of math operators to the package. The simplification function used by @oscarbenjamin is referenced via sympy.sqf:

@syms x::positive
fx = (x^6 + 2) / (8x^2)
dL = (1 + diff(fx, x)^2) |> sympy.sqf
2 * PI * integrate(fx * sqrt(dL), (x, 1, 3))

[rocco_sprmnt21]

it could also be done by hand

integrate(fx * sqrt(dL) = \frac{1}{16}(\frac{x^8}{8}+ \frac{3x^2}{2}+\frac{2x^{-4}}{-4})