Area of Surface of Revolution Integral too Hard to be Computed By JULIA’ SymPy and Python’ SymPy
68 views
Skip to first unread message
Freya the Goddess
Jan 15, 2023, 2:36:14 AM1/15/23
to sy...@googlegroups.com
Hi all,
I have a question: why SymPy (in JULIA and PYthon) unable to get the numerical answer for area of surface of revolution?
Is it impossible?
This is my question posted today on Julia Discourse:
--
С наилучшими пожеланиями, Богиня Фрейя
Atenciosamente, Freya the Goddess
Meilleurs voeux, Freya the Goddess
Liebe Grüße, Freya the Goddess
Best wishes, Freya the Goddess
よろしくお願いします、Freya the Goddess
最好的祝福，Freya the Goddess
Matakwa mema, Freya the Goddess
مع أطيب التمنيات ، فريا الإلهة
Aaron Meurer
Jan 16, 2023, 2:51:11 AM1/16/23
to sy...@googlegroups.com
If you want to compute a numerical answer using SymPy, use evalf() on the integral. But as others have pointed out if you are just doing a purely numerical integration you can also use the Julia functions to do that. SymPy will only be useful if you are interested in a symbolic answer, and in this case, a closed-form answer might not even exist.
Aaron Meurer
--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CALUh_%2B2dkMYNVgDvZHuZQW3%3DX9vS9EPqRv9fMVc5zbeUXYrmSw%40mail.gmail.com.
emanuel.c...@gmail.com
Jan 17, 2023, 4:51:20 AM1/17/23
to sympy
Note : in your question in Julia Discource, you state that your areao of revolution is 2*pi*integrate(f(x)*sqrt(1+f(x).diff(x)), (x, a, b)).
I beg to differ : it should be 2*pi*integrate(abs(f(x))*sqrt(1+abs(f(x).diff(x))), (x, a, b)). In your specific case, both f(x) and f(x).diff(x) are positive on [1 3], so there is no numerical difference.
BTW : Mathematica can get you an explicit expression for your integral, involving hypergeometric functions, but fails to compute its numerical values.
HTH,
HTH,
Oscar
Jan 21, 2023, 9:32:16 AM1/21/23
to sympy
On Sunday, 15 January 2023 at 07:36:14 UTC zaqhie...@gmail.com wrote:
Hi all,
I have a question: why SymPy (in JULIA and PYthon) unable to get the numerical answer for area of surface of revolution?
Is it impossible?
This is my question posted today on Julia Discourse:
Please ask questions here rather than posting a link to somewhere else.
You can numerically evaluate integrals using evalf:
In [1]: x = symbols("x")
...:
...: f = (x**6 + 2)/(8*x**2)
...: g = sqrt(1 + diff(f,x))
...:
...: h = 2*pi*Integral(((x**6 + 2)/(8*x**2))*sqrt(1 + diff(f,x)), (x, 1, 3))
In [2]: h
Out[2]:
3
⌠
⎮ ___________________
⎮ ╱ 3 6
⎮ ⎛ 6 ⎞ ╱ 3⋅x x + 2
⎮ ⎝x + 2⎠⋅ ╱ ──── + 1 - ──────
⎮ ╱ 4 3
⎮ ╲╱ 4⋅x
2⋅π⋅⎮ ────────────────────────────────── dx
⎮ 2
⎮ 8⋅x
⌡
1
In [3]: h.evalf()
Out[3]: 116.281297293490
--
Oscar
Kalevi Suominen
Jan 23, 2023, 12:22:01 PM1/23/23
to sympy
The derivative should actually be squared in the square root expression: sqrt(1 + f'(x)^2) (see e.g.
https://en.wikipedia.org/wiki/Surface_of_revolution), which then simplifies to a rational function (x^6 + 1)/(2*x^3) (unless I made a mistake).
Hence the integrand will be rational and SymPy should be able to handle it.
In general, the square root does simplify. In that case the result will be a hyperelliptic integral, which is non-elementary and cannot be
represented by means of common special functions. There is no support in SymPy for such integrals.
Kalevi Suominen
Kalevi Suominen
Jan 23, 2023, 12:23:45 PM1/23/23
to sympy
EDIT: In general, the square root does *not* simplify.
emanuel.c...@gmail.com
Jan 25, 2023, 4:14:59 PM1/25/23
to sympy
There seems to be some funamental differences between Sympy's integration algorithm an other ones (say Sage, Maxima, Giac, Fricas or Mathematica). Sympy :
```
>>> from sympy import *
>>> x=symbols("x", positive=True)
>>> f=Lambda((x), (x**6+2)/(8*x**2))
>>> integrate(sqrt(f(x).diff(x)**2+1)*f(x))
(Integral(2*Abs(x**4 - x**2 + 1)/x**5, x) + Integral(2*Abs(x**4 - x**2 + 1)/x**3, x) + Integral(x*Abs(x**4 - x**2 + 1), x) + Integral(x**3*Abs(x**4 - x**2 + 1), x))/16
>>> x=symbols("x", positive=True)
>>> f=Lambda((x), (x**6+2)/(8*x**2))
>>> integrate(sqrt(f(x).diff(x)**2+1)*f(x))
(Integral(2*Abs(x**4 - x**2 + 1)/x**5, x) + Integral(2*Abs(x**4 - x**2 + 1)/x**3, x) + Integral(x*Abs(x**4 - x**2 + 1), x) + Integral(x**3*Abs(x**4 - x**2 + 1), x))/16
```
Sage (either via Maxima's or Giac's integrators) :
```
sage: x=var("x", domain="positive")
sage: f(x)=(x^6+2)/(8*x^2)
sage: integrate((f(x).diff(x)^2+1).sqrt()*f(x),x)
1/128*x^8 + 1/32*x^2 + 1/32*(2*x^6 - 1)/x^4
```
sage: f(x)=(x^6+2)/(8*x^2)
sage: integrate((f(x).diff(x)^2+1).sqrt()*f(x),x)
1/128*x^8 + 1/32*x^2 + 1/32*(2*x^6 - 1)/x^4
```
BTW, calling these integrators from Sage :
```
sage: [[(u/v).full_simplify() for v in L] for u in L]
[[1, 1, 1, 1, (x^6 + 1)/sqrt(x^12 + 2*x^6 + 1)],
[1, 1, 1, 1, (x^6 + 1)/sqrt(x^12 + 2*x^6 + 1)],
[1, 1, 1, 1, (x^6 + 1)/sqrt(x^12 + 2*x^6 + 1)],
[1, 1, 1, 1, (x^6 + 1)/sqrt(x^12 + 2*x^6 + 1)],
[sqrt(x^12 + 2*x^6 + 1)/(x^6 + 1),
sqrt(x^12 + 2*x^6 + 1)/(x^6 + 1),
sqrt(x^12 + 2*x^6 + 1)/(x^6 + 1),
sqrt(x^12 + 2*x^6 + 1)/(x^6 + 1),
1]]
sage: S=integrate((f(x).diff(x)^2+1).sqrt()*f(x),x) ; S
1/128*x^8 + 1/32*x^2 + 1/32*(2*x^6 - 1)/x^4
sage: G=integrate((f(x).diff(x)^2+1).sqrt()*f(x),x, algorithm="giac") ; G
1/128*x^8*sgn(x^9 + x^3) + 3/32*x^2*sgn(x^9 + x^3) - 1/32*sgn(x^9 + x^3)/x^4
sage: F=integrate((f(x).diff(x)^2+1).sqrt()*f(x),x, algorithm="fricas") ; F
1/128*(x^12 + 12*x^6 - 4)/x^4
sage: Ma=integrate((f(x).diff(x)^2+1).sqrt()*f(x),x, algorithm="maxima") ; Ma
1/128*x^8 + 1/32*x^2 + 1/32*(2*x^6 - 1)/x^4
sage: M=mathematica.Refine(mathematica.Integrate((f(x).diff(x)^2+1).sqrt()*f(x),x).sage(), x>0).sage() ; Ma
1/128*x^8 + 1/32*x^2 + 1/32*(2*x^6 - 1)/x^4
sage: L=[S,G,F,Ma,M]
sage: [[(u/v).full_simplify() for v in L] for u in L]
[[1, 1, 1, 1, (x^6 + 1)/sqrt(x^12 + 2*x^6 + 1)],
[1, 1, 1, 1, (x^6 + 1)/sqrt(x^12 + 2*x^6 + 1)],
[1, 1, 1, 1, (x^6 + 1)/sqrt(x^12 + 2*x^6 + 1)],
[1, 1, 1, 1, (x^6 + 1)/sqrt(x^12 + 2*x^6 + 1)],
[sqrt(x^12 + 2*x^6 + 1)/(x^6 + 1),
sqrt(x^12 + 2*x^6 + 1)/(x^6 + 1),
sqrt(x^12 + 2*x^6 + 1)/(x^6 + 1),
sqrt(x^12 + 2*x^6 + 1)/(x^6 + 1),
1]]
[[1, 1, 1, 1, (x^6 + 1)/sqrt(x^12 + 2*x^6 + 1)],
[1, 1, 1, 1, (x^6 + 1)/sqrt(x^12 + 2*x^6 + 1)],
[1, 1, 1, 1, (x^6 + 1)/sqrt(x^12 + 2*x^6 + 1)],
[1, 1, 1, 1, (x^6 + 1)/sqrt(x^12 + 2*x^6 + 1)],
[sqrt(x^12 + 2*x^6 + 1)/(x^6 + 1),
sqrt(x^12 + 2*x^6 + 1)/(x^6 + 1),
sqrt(x^12 + 2*x^6 + 1)/(x^6 + 1),
sqrt(x^12 + 2*x^6 + 1)/(x^6 + 1),
1]]
sage: S=integrate((f(x).diff(x)^2+1).sqrt()*f(x),x) ; S
1/128*x^8 + 1/32*x^2 + 1/32*(2*x^6 - 1)/x^4
sage: G=integrate((f(x).diff(x)^2+1).sqrt()*f(x),x, algorithm="giac") ; G
1/128*x^8*sgn(x^9 + x^3) + 3/32*x^2*sgn(x^9 + x^3) - 1/32*sgn(x^9 + x^3)/x^4
sage: F=integrate((f(x).diff(x)^2+1).sqrt()*f(x),x, algorithm="fricas") ; F
1/128*(x^12 + 12*x^6 - 4)/x^4
sage: Ma=integrate((f(x).diff(x)^2+1).sqrt()*f(x),x, algorithm="maxima") ; Ma
1/128*x^8 + 1/32*x^2 + 1/32*(2*x^6 - 1)/x^4
sage: M=mathematica.Refine(mathematica.Integrate((f(x).diff(x)^2+1).sqrt()*f(x),x).sage(), x>0).sage() ; Ma
1/128*x^8 + 1/32*x^2 + 1/32*(2*x^6 - 1)/x^4
sage: L=[S,G,F,Ma,M]
sage: [[(u/v).full_simplify() for v in L] for u in L]
[[1, 1, 1, 1, (x^6 + 1)/sqrt(x^12 + 2*x^6 + 1)],
[1, 1, 1, 1, (x^6 + 1)/sqrt(x^12 + 2*x^6 + 1)],
[1, 1, 1, 1, (x^6 + 1)/sqrt(x^12 + 2*x^6 + 1)],
[1, 1, 1, 1, (x^6 + 1)/sqrt(x^12 + 2*x^6 + 1)],
[sqrt(x^12 + 2*x^6 + 1)/(x^6 + 1),
sqrt(x^12 + 2*x^6 + 1)/(x^6 + 1),
sqrt(x^12 + 2*x^6 + 1)/(x^6 + 1),
sqrt(x^12 + 2*x^6 + 1)/(x^6 + 1),
1]]
```
It *seems* that all these integrators, *except Sympy's*, use the fact that `x>0`. But this might also come from different forms of the integrand. Sympy :
```
>>> sqrt(f(x).diff(x)**2+1)*f(x)
(x**6 + 2)*sqrt((3*x**3/4 - (x**6 + 2)/(4*x**3))**2 + 1)/(8*x**2)
(x**6 + 2)*sqrt((3*x**3/4 - (x**6 + 2)/(4*x**3))**2 + 1)/(8*x**2)
```
Sage :
```
sage: sqrt(f(x).diff(x)^2+1)*f(x)
1/32*(x^6 + 2)*sqrt((3*x^3 - (x^6 + 2)/x^3)^2 + 16)/x^2
1/32*(x^6 + 2)*sqrt((3*x^3 - (x^6 + 2)/x^3)^2 + 16)/x^2
```
HTH,
Oscar Benjamin
Jan 25, 2023, 6:29:11 PM1/25/23
to sy...@googlegroups.com
On Wed, 25 Jan 2023 at 21:15, emanuel.c...@gmail.com
<emanuel.c...@gmail.com> wrote:
>
> There seems to be some funamental differences between Sympy's integration algorithm an other ones (say Sage, Maxima, Giac, Fricas or Mathematica). Sympy :
>
> ```
> >>> from sympy import *
> >>> x=symbols("x", positive=True)
> >>> f=Lambda((x), (x**6+2)/(8*x**2))
> >>> integrate(sqrt(f(x).diff(x)**2+1)*f(x))
> (Integral(2*Abs(x**4 - x**2 + 1)/x**5, x) + Integral(2*Abs(x**4 - x**2 + 1)/x**3, x) + Integral(x*Abs(x**4 - x**2 + 1), x) + Integral(x**3*Abs(x**4 - x**2 + 1), x))/16
> ```
>
> Sage (either via Maxima's or Giac's integrators) :
>
> ```
> sage: x=var("x", domain="positive")
> sage: f(x)=(x^6+2)/(8*x^2)
> sage: integrate((f(x).diff(x)^2+1).sqrt()*f(x),x)
> 1/128*x^8 + 1/32*x^2 + 1/32*(2*x^6 - 1)/x^4
Maybe it's because SymPy's integrate does not automatically factorise
<emanuel.c...@gmail.com> wrote:
>
> There seems to be some funamental differences between Sympy's integration algorithm an other ones (say Sage, Maxima, Giac, Fricas or Mathematica). Sympy :
>
> ```
> >>> from sympy import *
> >>> x=symbols("x", positive=True)
> >>> f=Lambda((x), (x**6+2)/(8*x**2))
> >>> integrate(sqrt(f(x).diff(x)**2+1)*f(x))
> (Integral(2*Abs(x**4 - x**2 + 1)/x**5, x) + Integral(2*Abs(x**4 - x**2 + 1)/x**3, x) + Integral(x*Abs(x**4 - x**2 + 1), x) + Integral(x**3*Abs(x**4 - x**2 + 1), x))/16
> ```
>
> Sage (either via Maxima's or Giac's integrators) :
>
> ```
> sage: x=var("x", domain="positive")
> sage: f(x)=(x^6+2)/(8*x^2)
> sage: integrate((f(x).diff(x)^2+1).sqrt()*f(x),x)
> 1/128*x^8 + 1/32*x^2 + 1/32*(2*x^6 - 1)/x^4
under the radical (I haven't checked the code to see if factorisation
is attempted):
In [29]: integrand = sqrt(f(x).diff(x)**2+1)*f(x)
In [30]: print(integrand)
(x**6 + 2)*sqrt((3*x**3/4 - (x**6 + 2)/(4*x**3))**2 + 1)/(8*x**2)
In [31]: print(integrand.factor())
(x**2 + 1)*(x**6 + 2)*Abs(x**4 - x**2 + 1)/(16*x**5)
In [32]: print(integrand.factor().integrate(x))
x**8/128 + 3*x**2/32 - 1/(32*x**4)
--
Oscar
Aaron Meurer
Jan 25, 2023, 8:04:31 PM1/25/23
to sy...@googlegroups.com
On Wed, Jan 25, 2023 at 4:29 PM Oscar Benjamin <oscar.j....@gmail.com> wrote:
On Wed, 25 Jan 2023 at 21:15, emanuel.c...@gmail.com
<emanuel.c...@gmail.com> wrote:
>
> There seems to be some funamental differences between Sympy's integration algorithm an other ones (say Sage, Maxima, Giac, Fricas or Mathematica). Sympy :
>
> ```
> >>> from sympy import *
> >>> x=symbols("x", positive=True)
> >>> f=Lambda((x), (x**6+2)/(8*x**2))
> >>> integrate(sqrt(f(x).diff(x)**2+1)*f(x))
> (Integral(2*Abs(x**4 - x**2 + 1)/x**5, x) + Integral(2*Abs(x**4 - x**2 + 1)/x**3, x) + Integral(x*Abs(x**4 - x**2 + 1), x) + Integral(x**3*Abs(x**4 - x**2 + 1), x))/16
> ```
>
> Sage (either via Maxima's or Giac's integrators) :
>
> ```
> sage: x=var("x", domain="positive")
> sage: f(x)=(x^6+2)/(8*x^2)
> sage: integrate((f(x).diff(x)^2+1).sqrt()*f(x),x)
> 1/128*x^8 + 1/32*x^2 + 1/32*(2*x^6 - 1)/x^4
Maybe it's because SymPy's integrate does not automatically factorise
under the radical (I haven't checked the code to see if factorisation
is attempted):
SymPy's integration algorithms that are currently implemented aren't particularly good with algebraic functions. They do work sometimes, but most of the time they don't. They are also quite sensitive to the exact form an expression is written in, and a typical algebraic function can be rewritten in multiple equivalent ways.
Aaron Meurer
In [29]: integrand = sqrt(f(x).diff(x)**2+1)*f(x)
In [30]: print(integrand)
(x**6 + 2)*sqrt((3*x**3/4 - (x**6 + 2)/(4*x**3))**2 + 1)/(8*x**2)
In [31]: print(integrand.factor())
(x**2 + 1)*(x**6 + 2)*Abs(x**4 - x**2 + 1)/(16*x**5)
In [32]: print(integrand.factor().integrate(x))
x**8/128 + 3*x**2/32 - 1/(32*x**4)
--
Oscar
--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxT5Yj1Wt%2BetUE-ep%3DrZqX%2BvZ5S-VQNyGas7ubL9kvY_uA%40mail.gmail.com.
Reply all
Reply to author
Forward
0 new messages