Sympy integration running into error, but Mathematica works fine

Question

I'm trying to integrate an equation using sympy, but the evaluation keeps erroring out with: TypeError: cannot add <class 'sympy.matrices.immutable.ImmutableDenseMatrix'> and <class 'sympy.core.numbers.Zero'>

However, when I integrate the same equation in Mathematica, I get the result I'm looking for. I'm hoping someone can explain what is happening under the hood with sympy and how its integration differs from mathematica's.

Because there are a lot of equations involved, I am only going to post the Python representation of the equations which are used to make up the variables called in the integration. To know for certain that the integration is an issue, I have independently checked each equation involved and made sure the results match between Python and Mathematica. I can confirm that only the integration fails.

Integration variable setup in Python:

import numpy as np
import sympy as sym
from sympy import I, Matrix, integrate
from sympy.integrals import Integral
from sympy.functions import exp
from sympy.physics.vector import cross, dot

c_speed = 299792458
lambda_u = 0.01
k_u = (2 * np.pi)/lambda_u
K_0 = 0.2
gamma_0 = 2.0
beta_0 = sym.sqrt(1 - 1 / gamma_0**2)
t = sym.symbols('t')
t_start = (-2 * lambda_u) / c_speed
t_end = (3 * lambda_u) / c_speed

beta_x_of_t = sym.Piecewise( (0.0, sym.Or(t < t_start, t > t_end)),
                             ((-1 * K_0)/gamma_0 * sym.sin(k_u * c_speed * t), sym.And(t_start <= t, t <= t_end)) )

beta_z_of_t = sym.Piecewise( (beta_0, sym.Or(t < t_start, t > t_end)),
                             (beta_0 * (1 - K_0**2/ (4 * beta_0 * gamma_0**2)) + K_0**2 / (4 * beta_0 * gamma_0**2) * sym.sin(2 * k_u * c_speed * t), sym.And(t_start <= t, t <= t_end)) )

beta_xp_of_t = sym.diff(beta_x_of_t, t)

beta_zp_of_t = sym.diff(beta_z_of_t, t)

Python Integration:

n = Matrix([(0, 0, 1)])
beta = Matrix([(beta_x_of_t, 0, beta_z_of_t)])
betap = Matrix([(beta_xp_of_t, 0, beta_zp_of_t)])

def rad(n, beta, betap):
  return integrate(n.cross((n-beta).cross(betap)), (t, t_start, t_end))

rad(n, beta, betap)
# Output is the error above

Mathematica integration:

rad[n_, beta_, betap_] :=
 NIntegrate[
  Cross[n, Cross[n - beta, betap]]
  , {t, tstart, tend}, AccuracyGoal -> 3]

rad[{0, 0, 1}, {betax[t], 0, betaz[t]}, {betaxp[t], 0, betazp[t]}]
# Output is {0.00150421, 0., 0.}

I did see similar question such as this one and this question, but I'm not entirely sure if they are relevant here (the second link seems to be a closer match, but not quite a fit).

Show some of your debugging efforts. I admit full tracebacks in sympy can be obscure, but that's a starting point for you - and us (if you want help!). If `integrate` returns an expression with an unevaluated integral, then you can legitimately say sympy can't integrate. With an error like yours, the problem probably lies in the argument expression. Which brings up another complaint - don't ask for sympy help without showing at least some of the intermediate expressions. I have no idea of what is being passed to `integrate`. I can't do that kind of math in my head! — hpaulj, Oct 25 '21 at 15:14
Thanks @hpaulj for your comment! This actually is my debugged example. The actual equation I am trying to integrate contains three additional parts, but I have narrowed down the error to this cross-product. I'll see if I can edit this to express some more. I'm trying to calculate radiation for one particle for a plasma — ml50, Oct 25 '21 at 19:28

JohanC · Accepted Answer · 2021-10-25T19:33:21.377

As you're working with imprecise floats (and even use np.pi instead of sym.pi), while sympy tries to find exact symbolic solutions, sympy's expressions get rather wild with constants like 6.67e-11 mixed with much larger values.

Here is an attempt to use more symbolic expressions, and only integrate the x-coordinate of the function.

import sympy as sym
from sympy import I, Matrix, integrate, S
from sympy.integrals import Integral
from sympy.functions import exp
from sympy.physics.vector import cross, dot

c_speed = 299792458
lambda_u = S(1) / 100
k_u = (2 * sym.pi) / lambda_u
K_0 = S(2) / 100
gamma_0 = 2
beta_0 = sym.sqrt(1 - S(1) / gamma_0 ** 2)
t = sym.symbols('t')
t_start = (-2 * lambda_u) / c_speed
t_end = (3 * lambda_u) / c_speed

beta_x_of_t = sym.Piecewise((0, sym.Or(t < t_start, t > t_end)),
                            ((-1 * K_0) / gamma_0 * sym.sin(k_u * c_speed * t), sym.And(t_start <= t, t <= t_end)))

beta_z_of_t = sym.Piecewise((beta_0, sym.Or(t < t_start, t > t_end)),
                            (beta_0 * (1 - K_0 ** 2 / (4 * beta_0 * gamma_0 ** 2)) + K_0 ** 2 / (
                                        4 * beta_0 * gamma_0 ** 2) * sym.sin(2 * k_u * c_speed * t),
                             sym.And(t_start <= t, t <= t_end)))

beta_xp_of_t = sym.diff(beta_x_of_t, t)

beta_zp_of_t = sym.diff(beta_z_of_t, t)

n = Matrix([(0, 0, 1)])
beta = Matrix([(beta_x_of_t, 0, beta_z_of_t)])
betap = Matrix([(beta_xp_of_t, 0, beta_zp_of_t)])

func = n.cross((n - beta).cross(betap))[0]

print(integrate(func, (t, t_start, t_end)))

This doesn't give an error, and outputs just zero.

Lambdify can be used to convert the function to numpy, and plot it.

func_np = sym.lambdify(t, func)

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

t_start_np = float(t_start.evalf())
t_end_np = float(t_end.evalf())

ts = np.linspace(t_start_np, t_end_np, 100000)
plt.plot(ts, func_np(ts))
print("approximate integral (np.trapz):", np.trapz(ts, func_np(ts)))
print("approximate integral (scipy's quad):", quad(func_np, t_start_np, t_end_np))

This outputs:

approximate integral (np.trapz): 0.04209721470548062
approximate integral (scipy's quad): (-2.3852447794681098e-18, 5.516333374450447e-10)

Note the huge values on the y-axis, and the small values on the x-axis. These values, as well as matematica's, could well be rounding errors.

Thanks, this really helps. As mentioned in my comment above, the original equation I am trying to integrate is much more complex, but breaking it up into the x, y, and z components seems to work for the y and z components. The x component seems to be taking a long time to integrate though (over 15 minutes and still running) so I might just switch to using scipy's quads function. Your answer definitely helped me understand how Sympy integrates though which was my goal. Marking as solved — ml50, Oct 25 '21 at 23:09

Sympy integration running into error, but Mathematica works fine

1 Answers1