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

3 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)

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

println("Area with QuadGK = ", dq)


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

1 Like

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.

2 Likes

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)
# 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())

"""



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


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})