[Python] Python Sympy vs Matlab Discrepancy in Calculating a Triple Integral

Can anyone think of a reason for the following major discrepancy (Python Sympy vs Matlab) in calculating the triple integral:



Both integrals should return the following result:



As seen below, Matlab symbolic toolbox does the calculation correctly:

syms l1 l2 l3 ell r theta phi real positive

% Integral A (Cartesian coordinates)
integrand_A = (1 - l1)*(1 - l2)*(1 - l3);
A = int(int(int(integrand_A, ...
l1, 0, sqrt(ell^2 - l2^2 - l3^2)), ...
l2, 0, sqrt(ell^2 - l3^2)), ...
l3, 0, ell);

% Integral B (Spherical coordinates)
x = r * sin(theta) * cos(phi);
y = r * sin(theta) * sin(phi);
z = r * cos(theta);
integrand_B = (1 - x)*(1 - y)*(1 - z) * r^2 * sin(theta);
B = int(int(int(integrand_B, ...
phi, 0, pi/2), ...
theta, 0, pi/2), ...
r, 0, ell);

% Expand and display results for comparison
A = expand(A);
B = expand(B);
disp('Result from Cartesian coordinates:')
pretty(A)
disp('Result from Spherical coordinates:')
pretty(B)
disp('Are they equal?')
isequal = logical(A == B);
disp(isequal)


Result from Cartesian coordinates:
6 5 4 3
ell ell pi ell 3 pi ell
- ---- + ---- - --------- + -------
48 5 16 6

Result from Spherical coordinates:
6 5 4 3
ell ell pi ell 3 pi ell
- ---- + ---- - --------- + -------
48 5 16 6

Are they equal?
1


However, Python-Sympy makes a serious mistake (B is correct but A is OFF):

import sympy as sp

# Define symbolic variables
l1, l2, l3, ell = sp.symbols('l1 l2 l3 ell', real=True, positive=True)
r, theta, phi = sp.symbols('r theta phi', real=True, positive=True)

# Integral A (Cartesian coordinates)
integrand_A = (1 - l1) * (1 - l2) * (1 - l3)
A = sp.integrate(integrand_A,
(l1, 0, sp.sqrt(ell**2 - l2**2 - l3**2)),
(l2, 0, sp.sqrt(ell**2 - l3**2)),
(l3, 0, ell))

# Integral B (Spherical coordinates)
x = r * sp.sin(theta) * sp.cos(phi)
y = r * sp.sin(theta) * sp.sin(phi)
z = r * sp.cos(theta)
integrand_B = (1 - x) * (1 - y) * (1 - z) * r**2 * sp.sin(theta)
B = sp.integrate(integrand_B,
(phi, 0, sp.pi/2),
(theta, 0, sp.pi/2),
(r, 0, ell))

# Display results
print("Result from Cartesian coordinates:")
print(A)
print("\nResult from Spherical coordinates:")
print(B)
print("\nAre they equal?")
print(A.equals(B))


Result from Cartesian coordinates:
-ell**6/48 + ell**5/5 + 3*pi*ell**4/16 - pi*ell**3/6

Result from Spherical coordinates:
-ell**6/48 + ell**5/5 - 3*pi*ell**4/16 + pi*ell**3/6

Are they equal?
False


Could there be a bug in Pythons' Sympy? I would appreciate any explanation or comment on the issue.

Continue reading...