Using eigenvalues to determine critical points with sympy

Question

I have this code which aims to solve a system of first order differential equations. It returns a symbolic solution.

import sympy as sym 
sym.init_printing()
from IPython.display import display_latex
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
def solve_sys(a):
    t=sym.symbols('t')
    x1=sym.Function('x1')
    y2=sym.Function('y2')
    u=sym.Function('u')
    v=sym.Function('v')
    eq1=sym.Eq(x1(t).diff(t),y2(t))
    eq2=sym.Eq(y2(t).diff(t),-x1(t)+a*(y2(t)-((y2(t)**3)/3)))
    matrix=sym.Matrix([eq1.rhs,eq2.rhs])
    matJ=matrix.jacobian([x1(t),y2(t)])
    lin_mat = matJ.subs({x1(t):0,y2(t):0})
    lin_mat*sym.Matrix([u(t),v(t)])
    evect=lin_mat.eigenvects()
    evals = list(lin_mat.eigenvals().keys())
    return evect, evals

I then have a function which aims to take the eigenvalues, perform some tests on them and then return the type of critical point they produce.

ef critical_point_test(a):
    blah,evals=solve_sys(a)
    if isinstance(evals[0],complex)==False and isinstance(evals[0],complex)==False and np.sign(evals[0])==np.sign(evals[1]) and evals[0]!=evals[1]:
        print('Node')
    elif isinstance(evals[0],complex)==False and isinstance(evals[0],complex)==False and np.sign(evals[0])!=np.sign(evals[1]):
        print('Saddle Point')
    elif evals[0]==evals[1]:
        print('Proper or Impropper Node')
    elif isinstance(evals[0],complex)==True and isinstance(evals[0],complex)==True and np.real(evals[0])!=0 and np.real(evals[1])!=0:
        print('Spiral Point')
    else:
        print('Center')

I get the error:

Invalid comparison of complex 1/2 - sqrt(3)*I/2

My question is how do I modify my solve_sys(a) function so that it returns a complex number that can then be interpreted by my second function. Thanks!

Answer 1

Using a = Symbol('a') and calculating the eigenvalues I get:

In [22]: evals                                                                                                                    
Out[22]: 
⎡      _______   _______        _______   _______⎤
⎢a   ╲╱ a - 2 ⋅╲╱ a + 2   a   ╲╱ a - 2 ⋅╲╱ a + 2 ⎥
⎢─ - ───────────────────, ─ + ───────────────────⎥
⎣2            2           2            2         ⎦

In [23]: e1, e2 = evals                                                                                                           

In [24]: e1                                                                                                                       
Out[24]: 
      _______   _______
a   ╲╱ a - 2 ⋅╲╱ a + 2 
─ - ───────────────────
2            2 

If you substitute for a then you can use is_real to find out if the eigenvalues are real:

In [25]: e1.subs(a, 1)                                                                                                            
Out[25]: 
1   √3⋅ⅈ
─ - ────
2    2  

In [26]: e1.subs(a, 1).is_real                                                                                                    
Out[26]: False

If it is real you can use is_positive to find out if it is positive:

In [27]: e1.subs(a, 4)                                                                                                            
Out[27]: 2 - √3

In [28]: e1.subs(a, 4).is_positive                                                                                                
Out[28]: True

Likewise is_zero tells you if the eigenvalue is zero.

Working symbolically without substituting you can ask for what values of a it will be positive:

In [29]: solve(e1 >= 0, a)                                                                                                         
Out[29]: 2 ≤ a

You can find for what value of a the eigenvalues will be exactly equal

In [39]: solve(Eq(e1, e2), a)                                                                                                     
Out[39]: [2]