使用特征值通过 sympy 确定临界点

Question

我有这段代码，旨在求解一阶微分方程组。它 return 是一个象征性的解决方案。

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

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

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) 函数，使其 return 成为一个复数，然后可以由我的第二个函数解释。谢谢！

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

如果用 a 代替，则可以使用 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]

使用特征值通过 sympy 确定临界点

Using eigenvalues to determine critical points with sympy

sympy

python-3.x

differential-equations