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

**URL:** <https://discourse.julialang.org/t/area-of-surface-of-revolution-integral-too-hard-to-be-computed-by-julia-sympy-and-python-sympy/92981>\
**Category:** General Usage\
**Tags:** question\
**Created:** [January 15, 2023, 5:24am UTC](https://discourse.julialang.org/t/area-of-surface-of-revolution-integral-too-hard-to-be-computed-by-julia-sympy-and-python-sympy/92981 "2023-01-15T05:24:05Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![Freya\_the\_Goddess](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/freya_the_goddess/32/36835_2.png) [@Freya\_the\_Goddess](https://discourse.julialang.org/u/Freya_the_Goddess)\
**Post date:** [January 15, 2023, 5:24am UTC](https://discourse.julialang.org/t/area-of-surface-of-revolution-integral-too-hard-to-be-computed-by-julia-sympy-and-python-sympy/92981/1 "2023-01-15T05:24:05Z")

</div>

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:

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/b/8/b8eeb0f5d40efebb803dcdcf024413bb6540904b.png)

with Python 3.9:

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/0/5/055d22ed3f325911b1748847c1ad17ca689f6370.png)

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

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [January 15, 2023, 5:30am UTC](https://discourse.julialang.org/t/area-of-surface-of-revolution-integral-too-hard-to-be-computed-by-julia-sympy-and-python-sympy/92981/2 "2023-01-15T05:30:57Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Ininterrompue](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ininterrompue/32/5594_2.png) [@Ininterrompue](https://discourse.julialang.org/u/Ininterrompue)\
**Post date:** [January 15, 2023, 5:42am UTC](https://discourse.julialang.org/t/area-of-surface-of-revolution-integral-too-hard-to-be-computed-by-julia-sympy-and-python-sympy/92981/3 "2023-01-15T05:42:56Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Freya\_the\_Goddess](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/freya_the_goddess/32/36835_2.png) [@Freya\_the\_Goddess](https://discourse.julialang.org/u/Freya_the_Goddess)\
**Post date:** [January 15, 2023, 5:54am UTC](https://discourse.julialang.org/t/area-of-surface-of-revolution-integral-too-hard-to-be-computed-by-julia-sympy-and-python-sympy/92981/4 "2023-01-15T05:54:27Z")

</div>

Thanks for the mailing list

---

<div class="post-metadata">

**Author:** ![uniment](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/uniment/32/24532_2.png) [@uniment](https://discourse.julialang.org/u/uniment)\
**Post date:** [January 15, 2023, 7:58am UTC](https://discourse.julialang.org/t/area-of-surface-of-revolution-integral-too-hard-to-be-computed-by-julia-sympy-and-python-sympy/92981/5 "2023-01-15T07:58:19Z")

</div>

Let’s try numerical integration.

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

```julia
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
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](https://en.wikipedia.org/wiki/Dual_number#Differentiation) to automatically differentiate the function.

```julia
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
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
julia> using QuadGK

julia> 2π*quadgk(g, 1, 3)[1]
116.2812972934902

```

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

---

<div class="post-metadata">

**Author:** ![Freya\_the\_Goddess](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/freya_the_goddess/32/36835_2.png) [@Freya\_the\_Goddess](https://discourse.julialang.org/u/Freya_the_Goddess)\
**Post date:** [January 15, 2023, 9:34am UTC](https://discourse.julialang.org/t/area-of-surface-of-revolution-integral-too-hard-to-be-computed-by-julia-sympy-and-python-sympy/92981/6 "2023-01-15T09:34:23Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Freya\_the\_Goddess](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/freya_the_goddess/32/36835_2.png) [@Freya\_the\_Goddess](https://discourse.julialang.org/u/Freya_the_Goddess)\
**Post date:** [January 15, 2023, 10:41am UTC](https://discourse.julialang.org/t/area-of-surface-of-revolution-integral-too-hard-to-be-computed-by-julia-sympy-and-python-sympy/92981/7 "2023-01-15T10:41:53Z")

</div>

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

```julia
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**

---

<div class="post-metadata">

**Author:** ![oscarbenjamin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscarbenjamin/32/22868_2.png) [@oscarbenjamin](https://discourse.julialang.org/u/oscarbenjamin)\
**Post date:** [May 19, 2023, 1:08pm UTC](https://discourse.julialang.org/t/area-of-surface-of-revolution-integral-too-hard-to-be-computed-by-julia-sympy-and-python-sympy/92981/8 "2023-05-19T13:08:13Z")

</div>

> [@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:

```julia
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:

```julia
>>> 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):

```julia
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.

---

<div class="post-metadata">

**Author:** ![Freya\_the\_Goddess](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/freya_the_goddess/32/36835_2.png) [@Freya\_the\_Goddess](https://discourse.julialang.org/u/Freya_the_Goddess)\
**Post date:** [May 21, 2023, 5:54am UTC](https://discourse.julialang.org/t/area-of-surface-of-revolution-integral-too-hard-to-be-computed-by-julia-sympy-and-python-sympy/92981/9 "2023-05-21T05:54:55Z")

</div>

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:

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

"""

```

 ![333](https://global.discourse-cdn.com/julialang/original/3X/3/d/3db65472fc3df5de592b87ce1fce01519c5317ef.png)

---

<div class="post-metadata">

**Author:** ![j\_verzani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_verzani/32/8551_2.png) [@j\_verzani](https://discourse.julialang.org/u/j_verzani)\
**Post date:** [May 21, 2023, 12:00pm UTC](https://discourse.julialang.org/t/area-of-surface-of-revolution-integral-too-hard-to-be-computed-by-julia-sympy-and-python-sympy/92981/10 "2023-05-21T12:00:36Z")

</div>

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`:

```julia
@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))

```

---

<div class="post-metadata">

**Author:** ![rocco\_sprmnt21](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rocco_sprmnt21/32/20127_2.png) [@rocco\_sprmnt21](https://discourse.julialang.org/u/rocco_sprmnt21)\
**Post date:** [May 21, 2023, 1:04pm UTC](https://discourse.julialang.org/t/area-of-surface-of-revolution-integral-too-hard-to-be-computed-by-julia-sympy-and-python-sympy/92981/11 "2023-05-21T13:04:46Z")

</div>

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