SymPy's nullspace struggling when Mathematica succeeds with ease

Question

I have been trying to find the nullspace of a symbolic matrix A in SymPy using the command A.nullspace. Now the computation does not finish (or takes longer than I waited for) for the matrix I entered.

The strange thing is that Mathematica does this computation within a fraction of a second. I'll give the specific example I am using as Mathematica input code (should run)

L = ({{0, -I a Sqrt[b] Conjugate[g] Conjugate[w], 0, I a Sqrt[b] g w, 
 2 g2 + 2 b g z Conjugate[g], 0, 0, 0, 
 2 g2}, {-I a Sqrt[b] g w, -g2 - b g z Conjugate[g] - 
  I (-dD + b d g Conjugate[g]), -I DD, 0, I a Sqrt[b] g w, 0, 0, 
 0, 0}, {0, -I DD, I dD - g2, 0, 0, I a Sqrt[b] g w, 0, 0, 
 0}, {I a Sqrt[b] Conjugate[g] Conjugate[w], 0, 
 0, -g2 - b g z Conjugate[g] - 
  I (dD - b d g Conjugate[g]), -I a Sqrt[b] Conjugate[
   g] Conjugate[w], 0, I DD, 0, 0}, {0, 
 I a Sqrt[b] Conjugate[g] Conjugate[w], 
 0, -I a Sqrt[b] g w, -2 g2 - 2 b g z Conjugate[g], -I DD, 0, 
 I DD, 0}, {0, 0, I a Sqrt[b] Conjugate[g] Conjugate[w], 
 0, -I DD, -2 g2 + I b d g Conjugate[g] - b g z Conjugate[g], 0, 
 0, I DD}, {0, 0, 0, I DD, 0, 
 0, -I dD - g2, -I a Sqrt[b] Conjugate[g] Conjugate[w], 0}, {0, 0,
  0, 0, I DD, 
 0, -I a Sqrt[b] g w, -2 g2 - I b d g Conjugate[g] - 
  b g z Conjugate[g], -I DD}, {0, 0, 0, 0, 0, I DD, 
 0, -I DD, -2 g2}});

MatrixForm[L]
NullSpace[L]

As said for the same matrix SymPy is struggling. Is there any reason for this? Is there any other way to solve this problem in python?

---

Here is my SymPy code:

import numpy as np
import scipy
import sympy as sp

a = sp.symbols("a", real=True, positive=True)
b = sp.symbols("b", real=True, positive=True)
g = sp.symbols("g", real=True, positive=True)
g2 = sp.symbols("g2", real=True, positive=True)
DD = sp.symbols("DD", real=True, positive=True)

dD = sp.symbols("dD", real=True)

w = sp.symbols("w")
z = sp.symbols("z")
d = sp.symbols("d")

L = sp.Matrix([[0, -1j*a*sp.sqrt(b)*sp.conjugate(g)*sp.conjugate(w), 0,            1j*a*sp.sqrt(b)*g*w, 
2*g2 + 2*b*g*z*sp.conjugate(g), 0, 0, 0, 
2*g2], [-1j*a*sp.sqrt(b)*g*w, -g2 - b*g*z*sp.conjugate(g) - 
1j*(-dD + b*d*g*sp.conjugate(g)), -1j*DD, 0, 1j*a*sp.sqrt(b)*g*w, 0, 0, 
0, 0], [0, -1j*DD, 1j*dD - g2, 0, 0, 1j*a*sp.sqrt(b)*g*w, 0, 0, 
0], [1j*a*sp.sqrt(b)*sp.conjugate(g)*sp.conjugate(w), 0, 
0, -g2 - b*g*z*sp.conjugate(g) - 
1j*(dD - b*d*g*sp.conjugate(g)), -1j*a*sp.sqrt(b)*sp.conjugate(
g)*sp.conjugate(w), 0, 1j*DD, 0, 0], [0, 
1j*a*sp.sqrt(b)*sp.conjugate(g)*sp.conjugate(w), 
0, -1j*a*sp.sqrt(b)*g*w, -2*g2 - 2*b*g*z*sp.conjugate(g), -1j*DD, 0, 
1j*DD, 0], [0, 0, 1j*a*sp.sqrt(b)*sp.conjugate(g)*sp.conjugate(w), 
0, -1j*DD, -2*g2 + 1j*b*d*g*sp.conjugate(g) - b*g*z*sp.conjugate(g), 0, 
0, 1j*DD], [0, 0, 0, 1j*DD, 0, 
0, -1j*dD - g2, -1j*a*sp.sqrt(b)*sp.conjugate(g)*sp.conjugate(w), 0], [0, 0,
0, 0, 1j*DD, 
0, -1j*a*sp.sqrt(b)*g*w, -2*g2 - 1j*b*d*g*sp.conjugate(g) - 
b*g*z*sp.conjugate(g), -1j*DD], [0, 0, 0, 0, 0, 1j*DD, 
0, -1j*DD, -2*g2]]);

sp.pprint(L)
sp.pprint(L.nullspace())

According to a comment, the following inforamtion may be relevant too:

• Real, positive parameters: a, b, g, g2, DD
• Real parameters: dD
• Complex parameters: w, z, d

Answer 1

The imaginary unit in sympy is sympy.I (not 1j). Apparently, 1j is acceptable by sympy, however, it is interpreted as a numercal quantity, which should be avoided in symbolic computations.

The following code gives the nullspace fairly quick. Note that I have also provided symbol attributes. These usually help sympy provide faster answers

a = sp.symbols("a", positive = True)
b = sp.symbols("b", positive = True)
g = sp.symbols("g", positive = True)
w = sp.symbols("w")
z = sp.symbols("z")
g2 = sp.symbols("g2", positive = True)
DD = sp.symbols("DD", positive = True)
dD = sp.symbols("dD", real = True)
d = sp.symbols("d")

L = sp.Matrix([[0, -sp.I*a*sp.sqrt(b)*sp.conjugate(g)*sp.conjugate(w), 0,            sp.I*a*sp.sqrt(b)*g*w, 
2*g2 + 2*b*g*z*sp.conjugate(g), 0, 0, 0, 
2*g2], [-sp.I*a*sp.sqrt(b)*g*w, -g2 - b*g*z*sp.conjugate(g) - 
sp.I*(-dD + b*d*g*sp.conjugate(g)), -sp.I*DD, 0, sp.I*a*sp.sqrt(b)*g*w, 0, 0, 
0, 0], [0, -sp.I*DD, sp.I*dD - g2, 0, 0, sp.I*a*sp.sqrt(b)*g*w, 0, 0, 
0], [sp.I*a*sp.sqrt(b)*sp.conjugate(g)*sp.conjugate(w), 0, 
0, -g2 - b*g*z*sp.conjugate(g) - 
sp.I*(dD - b*d*g*sp.conjugate(g)), -sp.I*a*sp.sqrt(b)*sp.conjugate(
g)*sp.conjugate(w), 0, sp.I*DD, 0, 0], [0, 
sp.I*a*sp.sqrt(b)*sp.conjugate(g)*sp.conjugate(w), 
0, -sp.I*a*sp.sqrt(b)*g*w, -2*g2 - 2*b*g*z*sp.conjugate(g), -sp.I*DD, 0, 
sp.I*DD, 0], [0, 0, sp.I*a*sp.sqrt(b)*sp.conjugate(g)*sp.conjugate(w), 
0, -sp.I*DD, -2*g2 + sp.I*b*d*g*sp.conjugate(g) - b*g*z*sp.conjugate(g), 0, 
0, sp.I*DD], [0, 0, 0, sp.I*DD, 0, 
0, -sp.I*dD - g2, -sp.I*a*sp.sqrt(b)*sp.conjugate(g)*sp.conjugate(w), 0], [0, 0,
0, 0, sp.I*DD, 
0, -sp.I*a*sp.sqrt(b)*g*w, -2*g2 - sp.I*b*d*g*sp.conjugate(g) - 
b*g*z*sp.conjugate(g), -sp.I*DD], [0, 0, 0, 0, 0, sp.I*DD, 
0, -sp.I*DD, -2*g2]]);

Ln = L.nullspace()