zero in sympy: Issue #22425

[Zoufiné Lauer-Baré]

Dear all,

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)...

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]

On Thu, 4 Nov 2021 at 23:13, Zoufiné Lauer-Baré <zoufin...@gmail.com> wrote:
>
> 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)...

I'm not sure what the issue here is:

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é]

Hello Oscar,

 thank you for the fast answer. I'll test the chop option. Further, I commented on this idea in "Zero in SymPy #22425".

Regards,

Zoufiné

--
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.