0

我正在尝试解决 sympy 中 Sturm Liouville 问题的零点

这是我的代码

from sympy import *                                                            
init_printing(use_latex='mathjax')                                             

# Function Definitions                                                         
u = Function('u')                                                              
X = Function('X')                                                              

# Variables and Dummies                                                        
x  = Dummy('x')                                                                                                                                               
mu  = Symbol('mu', nonzero=True, positive=True)                                
C1 = Symbol('C1', nonzero=True)                                                
C2 = Symbol('C2', nonzero=True)                                                

print("STURM LIOUVILLE PROBLEM IN X")                                        
print("GENERAL SOLUTION")                                                      
# --- --- ---#                                                                 
print("")                                                                      
gen_sol = dsolve(Eq(Derivative(X(x), x, 2) + (mu**2)*X(x), 0), X(x))           
print(pretty(gen_sol))                                                         
print("")                                                                      
print("SOLVE FOR BOUNDARY")                                                    
bc_2 = gen_sol.rhs.diff(x).subs(x,5).subs(Symbol('C1')*mu, Symbol('C2'))       
print("DETERMINE THAT C2 != 0 (NON TRIVIAL)")                                  
print(pretty(bc_2))                                                            
print("")                                                                      
print("SOLVE FOR REMAINING SOLUTION")                                          
bc_2 = bc_2.subs(Symbol('C2'), 1)                                              
print(pretty(bc_2))                                                            

print("")                                                                      
linear_solutions = solve([bc_2], [mu])                                         
print(pretty(linear_solutions))                               

不幸的是,我收到以下 NotImplemented 错误。

有没有另一种方法可以解决可以处理此操作的 sympy 中的零?

STURM LIOUVILLE PROBLEM IN X
GENERAL SOLUTION

X(x) = C₁⋅sin(x⋅μ) + C₂⋅cos(x⋅μ)

SOLVE FOR BOUNDARY
DETERMINE THAT C2 != 0 (NON TRIVIAL)
-C₂⋅μ⋅sin(5⋅μ) + C₂⋅cos(5⋅μ)

SOLVE FOR REMAINING SOLUTION
-μ⋅sin(5⋅μ) + cos(5⋅μ)

Traceback (most recent call last):
  File "stack_overflow.py", line 31, in <module>
    linear_solutions = solve([bc_2], [mu])
  File "/home/cgould/.local/lib/python3.8/site-packages/sympy/solvers/solvers.py", line 1176, in solve
    solution = _solve_system(f, symbols, **flags)
  File "/home/cgould/.local/lib/python3.8/site-packages/sympy/solvers/solvers.py", line 1943, in _solve_system
    raise NotImplementedError('could not solve %s' % eq2)

(为清晰而编辑)这里的想法是我试图找到下面的零

-μ⋅sin(5⋅μ) + cos(5⋅μ) = 0

大致翻译为:

cot(5μ) = μ 
4

1 回答 1

1

solve将只求解一个已经编程了封闭形式解的方程,给出一个精确的符号解。nsolve将找到一个数值近似值:

>>> μ = mu
>>> nsolve(-μ*sin(5*μ) + cos(5*μ), 0)
0.262767543298580

由于您的方程有多个根,因此您必须使用不同的初始猜测来获得不同的解(或使用二分法并给出您想要找到根的区域)。一个很好的方法是找到前两个根,然后使用它们之间的差异来猜测到下一个根的距离:一阶延续方法。

>>> r = []
>>> r.append(nsolve(-μ*sin(5*μ) + cos(5*μ),0))
>>> r.append(nsolve(-μ*sin(5*μ) + cos(5*μ),.5))
>>> for i in range(5):
...   r.append(nsolve(-μ*sin(5*μ) + cos(5*μ),2*r[-1] - r[-2]))
...
>>> [i.round(2) for i in r]
[0.26, 0.81, 1.38, 1.98, 2.59, 3.20, 3.82]
于 2020-04-27T04:35:33.973 回答