Sympy fails, wxMaxima not

Question

I am trying to solve the following indefinite integral both with wxMaxima and sympy:

integrate(r^2*sqrt(R^2-r^2),r)

In Maxima I do got an answer but not in sympy. I do not understand why. I am a heavy user of Python and will love to do symbolic math in Python, but since sympy did not solve this, I am still stuck with Maxima.

Is it something I am doing wrong or is Maxiam better? (I also solved the same in Mathematica)

I got the following answer in wxMaxima:

f:r^2*sqrt(R^2-r^2);
g:integrate(f,r);

gives this answer:

g:(R^4*asin(r/abs(R)))/8-(r*(R^2-r^2)^(3/2))/4+(r*R^2*sqrt(R^2-r^2))/8  

It is looking ugly, but forget that. The point here is that sympy cannot solve this integral. Trying to solve the same in sympy with this code:

import sympy as sy
import math
R,r = sy.symbols('R r')
g = sy.integrate(r**2*(R**2-r**2)**0.5,r)
print g

Gives the this error message:

Traceback (most recent call last):
  File "E:\UD\Software\BendStiffener\calculate-moment\moment-calculation-equations\sympy-test.py", line 4, in <module>
    g = sy.integrate(r**2*(R**2-r**2)**0.5,r)
  File "C:\Python27\lib\site-packages\sympy\utilities\decorator.py", line 35, in threaded_func
    return func(expr, *args, **kwargs)
  File "C:\Python27\lib\site-packages\sympy\integrals\integrals.py", line 1232, in integrate
    risch=risch, manual=manual)
  File "C:\Python27\lib\site-packages\sympy\integrals\integrals.py", line 487, in doit
    conds=conds)
  File "C:\Python27\lib\site-packages\sympy\integrals\integrals.py", line 876, in _eval_integral
    h = meijerint_indefinite(g, x)
  File "C:\Python27\lib\site-packages\sympy\integrals\meijerint.py", line 1596, in meijerint_indefinite
    res = _meijerint_indefinite_1(f.subs(x, x + a), x)
  File "C:\Python27\lib\site-packages\sympy\integrals\meijerint.py", line 1646, in _meijerint_indefinite_1
    r = hyperexpand(r.subs(t, a*x**b))
  File "C:\Python27\lib\site-packages\sympy\simplify\hyperexpand.py", line 2482, in hyperexpand
    return f.replace(hyper, do_replace).replace(meijerg, do_meijer)
  File "C:\Python27\lib\site-packages\sympy\core\basic.py", line 1351, in replace
    rv = bottom_up(self, rec_replace, atoms=True)
  File "C:\Python27\lib\site-packages\sympy\simplify\simplify.py", line 4082, in bottom_up
    rv = F(rv)
  File "C:\Python27\lib\site-packages\sympy\core\basic.py", line 1336, in rec_replace
    new = _value(expr, result)
  File "C:\Python27\lib\site-packages\sympy\core\basic.py", line 1280, in <lambda>
    _value = lambda expr, result: value(*expr.args)
  File "C:\Python27\lib\site-packages\sympy\simplify\hyperexpand.py", line 2479, in do_meijer
    allow_hyper, rewrite=rewrite)
  File "C:\Python27\lib\site-packages\sympy\simplify\hyperexpand.py", line 2375, in _meijergexpand
    t, 1/z0)
  File "C:\Python27\lib\site-packages\sympy\simplify\hyperexpand.py", line 2335, in do_slater
    resid = residue(integrand, s, b_ + r)
  File "C:\Python27\lib\site-packages\sympy\series\residues.py", line 57, in residue
    s = expr.series(x, n=0)
  File "C:\Python27\lib\site-packages\sympy\core\expr.py", line 2435, in series
    rv = self.subs(x, xpos).series(xpos, x0, n, dir, logx=logx)
  File "C:\Python27\lib\site-packages\sympy\core\expr.py", line 2442, in series
    s1 = self._eval_nseries(x, n=n, logx=logx)
  File "C:\Python27\lib\site-packages\sympy\core\mul.py", line 1446, in _eval_nseries
    terms = [t.nseries(x, n=n, logx=logx) for t in self.args]
  File "C:\Python27\lib\site-packages\sympy\core\expr.py", line 2639, in nseries
    return self._eval_nseries(x, n=n, logx=logx)
  File "C:\Python27\lib\site-packages\sympy\functions\special\gamma_functions.py", line 168, in _eval_nseries
    return super(gamma, self)._eval_nseries(x, n, logx)
  File "C:\Python27\lib\site-packages\sympy\core\function.py", line 598, in _eval_nseries
    term = e.subs(x, S.Zero)
  File "C:\Python27\lib\site-packages\sympy\core\basic.py", line 892, in subs
    rv = rv._subs(old, new, **kwargs)
  File "C:\Python27\lib\site-packages\sympy\core\basic.py", line 1006, in _subs
    rv = fallback(self, old, new)
  File "C:\Python27\lib\site-packages\sympy\core\basic.py", line 983, in fallback
    rv = self.func(*args)
  File "C:\Python27\lib\site-packages\sympy\core\function.py", line 382, in __new__
    return result.evalf(mlib.libmpf.prec_to_dps(pr))
  File "C:\Python27\lib\site-packages\sympy\core\evalf.py", line 1317, in evalf
    result = evalf(self, prec + 4, options)
  File "C:\Python27\lib\site-packages\sympy\core\evalf.py", line 1217, in evalf
    re, im = x._eval_evalf(prec).as_real_imag()
  File "C:\Python27\lib\site-packages\sympy\core\function.py", line 486, in _eval_evalf
    v = func(*args)
  File "C:\Python27\lib\site-packages\sympy\mpmath\ctx_mp_python.py", line 993, in f
    return ctx.make_mpf(mpf_f(x._mpf_, prec, rounding))
  File "C:\Python27\lib\site-packages\sympy\mpmath\libmp\gammazeta.py", line 1947, in mpf_gamma
    raise ValueError("gamma function pole")
ValueError: gamma function pole

Answer 1

In general, SymPy will perform better with rational numbers than with floating point numbers, especially when those numbers are powers.

This is because "nice" closed form solutions often only exist for exact powers. Consider for instance the difference between this

In [39]: integrate(r**2*sqrt(R**2-r**2), r)
Out[39]:
⎧     4      ⎛r⎞
⎪  ⅈ⋅R ⋅acosh⎜─⎟            3                      3                  5             │ 2│
⎪            ⎝R⎠         ⅈ⋅R ⋅r             3⋅ⅈ⋅R⋅r                ⅈ⋅r              │r │
⎪- ───────────── + ───────────────── - ───────────────── + ───────────────────  for │──│ > 1
⎪        8                 _________           _________             _________      │ 2│
⎪                         ╱       2           ╱       2             ╱       2       │R │
⎪                        ╱       r           ╱       r             ╱       r
⎪                  8⋅   ╱   -1 + ──    8⋅   ╱   -1 + ──    4⋅R⋅   ╱   -1 + ──
⎪                      ╱          2        ╱          2          ╱          2
⎪                    ╲╱          R       ╲╱          R         ╲╱          R
⎨
⎪     4     ⎛r⎞
⎪    R ⋅asin⎜─⎟          3                     3                 5
⎪           ⎝R⎠         R ⋅r              3⋅R⋅r                 r
⎪    ────────── - ──────────────── + ──────────────── - ──────────────────       otherwise
⎪        8                ________           ________             ________
⎪                        ╱      2           ╱      2             ╱      2
⎪                       ╱      r           ╱      r             ╱      r
⎪                 8⋅   ╱   1 - ──    8⋅   ╱   1 - ──    4⋅R⋅   ╱   1 - ──
⎪                     ╱         2        ╱         2          ╱         2
⎩                   ╲╱         R       ╲╱         R         ╲╱         R

and this

In [40]: integrate(r**2*(R**2-r**2)**0.5001, r)
Out[40]:
                                  ⎛             │  2  2⋅ⅈ⋅π⎞
                   1.0002  3  ┌─  ⎜-0.5001, 3/2 │ r ⋅ℯ     ⎟
0.333333333333333⋅R      ⋅r ⋅ ├─  ⎜             │ ─────────⎟
                             2╵ 1 ⎜    5/2      │      2   ⎟
                                  ⎝             │     R    ⎠ 

The power is almost 0.5, but the answer requires the use of a special function to represent. SymPy may not notice that 0.5 is supposed to be 1/2 (and the inexactness of floating point numbers doesn't help with this).

With that being said, I would consider this to be a SymPy bug, especially since it can compute integrate(r**2*(R**2-r**2)**0.5001, r).

Answer 2

You get a solution when you slightly rewrite your equation:

import sympy as sy
import math
R, r = sy.symbols('R r')

g = sy.integrate(r**2 * sy.sqrt((R**2 - r**2)), r)

print g.simplify()

So instead of using expr**0.5, I use sy.sqrt(expr). That gives

Piecewise((I*R*(-R**3*acosh(r/R) - R**2*r*sqrt((-R**2 + r**2)/R**2) + 2*r**3*sqrt((-R**2 + r**2)/R**2))/8, Abs(r**2/R**2) > 1), (R*(R**3*asin(r/R) - R**2*r*sqrt(1 - r**2/R**2) + 2*r**3*sqrt(1 - r**2/R**2))/8, True))

Whether that is the same as the results from Maxima is difficult to tell as sympy gives you the solution in two parts which depend on whether the argumnt in sqrt is larger or smaller than 0. You can try to use actual bounds and check whether you get the same results for the Maxima solution and the second part of the result sympy gives you.

You can access the second part of the solution like this:

g.args[1][0]

which gives you:

 R**4*asin(r/R)/8 - R**3*r/(8*sqrt(1 - r**2/R**2)) + 3*R*r**3/(8*sqrt(1 - r**2/R**2)) - r**5/(4*R*sqrt(1 - r**2/R**2))

You can also get the simplified version by doing:

 g.args[1][0].simplify()

which gives you:

R*(R**3*asin(r/R) - R**2*r*sqrt(1 - r**2/R**2) + 2*r**3*sqrt(1 - r**2/R**2))/8

That looks now quite similar to the result obtained by Maxima.