python-3.x - 使用特征值确定临界点与 sympy

Question

我有这段代码，旨在解决一阶微分方程组。它返回一个符号解决方案。

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

然后我有一个函数，旨在获取特征值，对它们执行一些测试，然后返回它们产生的临界点的类型。

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

我得到错误：

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

我的问题是如何修改我的 solve_sys(a) 函数，以便它返回一个复数，然后我的第二个函数可以解释该复数。谢谢！

标签： python-3.xsympydifferential-equations

Answer 1

使用a = Symbol('a')和计算我得到的特征值：

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 

如果您替换athen 您可以使用is_real来确定特征值是否为实数：

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

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

如果它是真实的，您可以使用is_positive它来确定它是否是正面的：

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

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

同样is_zero告诉您特征值是否为零。

在不替代的情况下象征性地工作，您可以询问a它的哪些值是积极的：

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

您可以找到特征值a的值将完全相等

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