zero in sympy: Issue #22425
Zoufiné Lauer-Baré
Sometimes SymPy hesitates to return zero... I've encountered this problem in three applications. There may be a solution to this already, however I haven't seen it yet.
- Problem 1: Real symmetric Matrices have only real eigenvalues...
- Problem 2: Analiticity of Möbius transform
- Problem 3: Stationary Points of Himmelblau Function
Problem 1:
A = sym.Matrix(([1, 4, -2],
[4, 0, 0],
[-2, 0, 3]))
should have only real eigenvalues, since it is symmetric, but SymPy returns complex eigenvalues with an imaginary part of the orrder 10**(-126)...
Problem 1:
The Moebius transform
f = (7.6sym.I(x + sym.Iy) - csym.I)/(-dx - dsym.I*y + 1)
fulfills the Cauchy-Rieman conditions, i.e.
sym.simplify(sym.diff(sym.im(f), y)-sym.diff(sym.re(f), x)) = 0
and
sym.simplify(sym.diff(sym.im(f), x)+sym.diff(sym.re(f), y)) = 0
However, when using numerical values for d and c
f = (7.6sym.I(x + sym.Iy) - 15.3215831575369sym.I)/(-2.01599778388644x - 2.01599778388644sym.I*y + 1)
it does not fullfill the Cauchy-Riemann conidtions anymore
Problem 3:
The Himmelblau function
f = (x**2+y-11)2+(x+y2-7)**2
has four stationary points in R x R (https://en.wikipedia.org/wiki/Himmelblau%27s_function)
However
system = [sym.diff(f, x),
sym.diff(f, y),
]
solSet = sym.nonlinsolve(system,[x,y])
solSetReal=[]
for i in list(solSet):
if i[0].is_real and i[1].is_real:
solSetReal.append(i)
solSetReal
returns only one stationary point in R x R.
While the imaginray parts of the other stationary points are actually zero.
sym.im(list(solSet)[1][0]).evalf()
gives a value of the order 10**(-125)...
I attached a Jupyter Notebbok with the three examples
Thank you!
Regards,
Zoufiné
Oscar Benjamin
>
> Dear all,
Hi and thanks for reporting this.
> I just created an issue in the SymPy github project (Zero in SymPy
> #22425 opened 2 minutes ago by zolabar). Here goes the content, may be someone knows this topic and there are already issues on this.
>
> Sometimes SymPy hesitates to return zero... I've encountered this problem in three applications. There may be a solution to this already, however I haven't seen it yet.
>
> Problem 1: Real symmetric Matrices have only real eigenvalues...
> Problem 2: Analiticity of Möbius transform
> Problem 3: Stationary Points of Himmelblau Function
>
> Problem 1:
>
> A = sym.Matrix(([1, 4, -2],
> [4, 0, 0],
> [-2, 0, 3]))
>
> should have only real eigenvalues, since it is symmetric, but SymPy returns complex eigenvalues with an imaginary part of the orrder 10**(-126)...
In [42]: A = sym.Matrix(([1, 4, -2],
...: [4, 0, 0],
...: [-2, 0, 3]))
In [43]: nroots(A.charpoly())
Out[43]: [-3.79943573866291, 2.29524145208425, 5.50419428657866]
In [44]: {e.evalf() for e in A.eigenvals()}
Out[44]: {-3.79943573866291 + 0.e-20⋅ⅈ, 2.29524145208425 + 0.e-20⋅ⅈ,
5.50419428657866 + 0.e-20⋅ⅈ}
The complex parts here show as zero. You can get rid of them
completely with chop:
In [45]: {e.evalf(chop=True) for e in A.eigenvals()}
Out[45]: {-3.79943573866291, 2.29524145208425, 5.50419428657866}
The source of the complex part is casus irreducibilis:
https://en.wikipedia.org/wiki/Casus_irreducibilis
That shows why it is better to compute numerical roots directly from
the polynomial rather than from a radical expression for the roots.
I don't immediately have time to investigate the other two points
(copy-pasting the code didn't work).
--
Oscar
Zoufiné Lauer-Baré
--
You received this message because you are subscribed to a topic in the Google Groups "sympy" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sympy/rMAWl_g8z9s/unsubscribe.
To unsubscribe from this group and all its topics, send an email to sympy+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxSTY4mUUV5cYVvRc77yKUURGKR%3DCOoBm0ZAhUtgxXkV%2Bg%40mail.gmail.com.