Symbolic solution of linear equation system returns wrong results

[hpaantee]

Related to #442. Solving a linear system Ax=b symbolically via A.LUsolve(b) returns a wrong result. Both when using .subs or Lambdify. The result contains values of up to e+100. When solving the linear system numerically I get the correct result, with no value exceeding e-02.

What causes this difference? In this comment (sympy/sympy#25208 (comment)) one contributor of sympy also mentions the accuracy difference. Where does this come from?

[hpaantee]

Small update: By converting every number in A and b from python's float to sympy's (and later symengine's) Float and specifying a precision via the dps argument, I got quite some improvement in the result. However, the result is still quite off and not yet usable. Even setting dps=1000 does not help.

[isuruf]

With symbolic values, you can't do pivoting correctly. LU decomposition using Gaussian elimination can be numerically unstable when small values are chosen as pivots and roundoff errors become larger.