1

I find a serious but unexplainable bug in the evalf function when evaluating a symbolic determinant of a matrix, see MWE below

import sympy
z1,z2,z3,mv1,mv2,mv3,R,T,D12,D23,D13,D14,D24,D34,x1,x2,x3 = sympy.symbols('z1,z2,z3,mv1,mv2,mv3,R,T,D12,D23,D13,D14,D24,D34,x1,x2,x3')
matrix = sympy.Matrix([ [  x2/D12 +  x3/D13 + 1/D14 ,         -x2/D12          ,           -x3/D13        ,  mv1/R/T  ] ,
                      [             -x1/D12         ,  x1/D12 + x3/D23 + 1/D24 ,           -x3/D23        ,  mv2/R/T  ] ,
                      [             -x1/D13         ,         -x2/D23          ,  x1/D13 + x2/D23 + 1/D34 ,  mv3/R/T  ] ,
                      [               z1*x1         ,           z2*x2          ,            z3*x3         ,     0     ]  ])

det1 = sympy.det(matrix.evalf(subs={ z1:-1,z2:-1,z3:1,mv1:50E-6,mv2:90E-6,mv3:150E-6,R:8.314,T:298,D12:1E-9,D23:9E-9,D13:3E-9,D14:0.165171466666667,D24:0.0917619259259259,D34:0.0550571555555556,x1:0.3,x2:0.2,x3:0.5}))
det2 = sympy.det(matrix).evalf(subs={z1:-1,z2:-1,z3:1,mv1:50E-6,mv2:90E-6,mv3:150E-6,R:8.314,T:298,D12:1E-9,D23:9E-9,D13:3E-9,D14:0.165171466666667,D24:0.0917619259259259,D34:0.0550571555555556,x1:0.3,x2:0.2,x3:0.5})
det3 = sympy.det(matrix).subs({      z1:-1,z2:-1,z3:1,mv1:50E-6,mv2:90E-6,mv3:150E-6,R:8.314,T:298,D12:1E-9,D23:9E-9,D13:3E-9,D14:0.165171466666667,D24:0.0917619259259259,D34:0.0550571555555556,x1:0.3,x2:0.2,x3:0.5})
print(det1)
print(det2)
print(det3)

produces answers of:

det1: 3.05615930451006e-7

det2: -439498375.118201

det3: 0

showing that evalf() on the symbolic output from sympy.det() has a serious bug -- all three answers should be near zero (to within ~1E-7 of numerical error). Method 1 and 2 work fine but clearly the det2 method has a huge bug. Using higher precision with evalf doesn't help.

Grateful for anyone to help: does anyone know if I've made a mistake here, or does sympy have a big bug here?

user3196617
  • 156
  • 2
  • 6
  • When I ran your code, I could not get `det2` to finish in a reasonable amount of time to check what value it gave (how long does it take to evaluate?). It seemed like it was the determinant rather than the evaluation slowing it down, so I used [this approach](https://stackoverflow.com/questions/37026935/speeding-up-computation-of-symbolic-determinant-in-sympy) to speed up the evaluation and once I did I got a value of `-6.25389203512624e-15` for `det2`. – Tyberius Sep 02 '21 at 15:23
  • That "det2 = ...." line requires ~3 seconds on my 2018 macbook pro. – user3196617 Sep 02 '21 at 15:46
  • I narrowed down the origin of the bug. The det2 line calculates better if run as: "sympy.det(matrix).subs({ z1:-1,z2:-1,z3:1}).evalf(subs={mv1:50E-6,mv2.........." Thus, evalf() has a bug handling the "z1:-1,z2:-1,z3:1". But what exactly is it's problem? – user3196617 Sep 02 '21 at 15:54
  • Weird, I'm on a pretty new Dell computer using Python3.8 from Anaconda. `det1` is done instantly, I have yet to get `det2` to finish after minutes of waiting. I tried updating sympy and that didn't seem to help. – Tyberius Sep 02 '21 at 15:59
  • If I change `det2` and `det3` to use `matrix.det(method="det_LU")` rather `sympy.det(matrix)`, they run much faster for me (basically instant) and I get `det2=-6.25389203512624e-15` and `det3=-4.76837158203125e-7`. If I increase the precision of `det1` and `det2` to 30, I get the same result for each of them to within a very high precision: `-6.253892e-15` – Tyberius Sep 02 '21 at 16:17
  • Ok I was using python 3.7 and sympy 1.7 or 1.6 (some version before 1.8). Just now I upgraded to sympy 1.8 and python 3.9, then tested the same code above from the orignal post, and yes I agree with you now: the det2() calculation freezes sympy indefinitely. – user3196617 Sep 02 '21 at 18:06
  • So it seems there was an old bug in sympy 1.7 (or whichever old version of sympy I was using) with the evalf() function, and now in sympy 1.8 there's a totally new and different bug with the evalf() function. – user3196617 Sep 02 '21 at 18:08

0 Answers0