python - 为 Sympy 解决方案建议解决方案

Question

我正在尝试使用 SymPy 的 solve 求解方程对于许多不同的输入值。这是到目前为止我的代码:

E = np.linspace(Emin,Emax,((Emax-Emin)/Emesh)+1)
z1 = np.zeros(len(E))+0j

for i in xrange(len(E)):    
    z = Symbol('z')
    z1[i] = solve( (d - E[i])/s - z - 1/z, z)

这里d和s是常数并且 E是一个一维数组。该代码工作正常，但速度非常慢。因为E[i]变化缓慢，z1[i]几乎等于z1[i+1] 。而不是每次为每个新的E[i]从头开始求解方程，我想使用结果 z1[i]作为解决 z1[i+1] 的起点。另一种选择是对 z[i+1] 的允许值进行限制。基于 z[i] 的值。基本上我希望尽可能加快这部分的速度。

有什么建议吗？

Answer 1

如果你看看你的等式，

(d - e) / s - z - 1/z == 0

乘以-z。

z**2 + ((e - d) / 2) * z + 1 == 0

现在这是一个我们可以轻松求解的二次方程，

b = (e - d) / s

z = (-b +/- (b*b - 4)**0.5) / 2

numpy 代码可以做到这一点，

import numpy as np

# some made-up numbers
e_min = 1.9
e_max = 2.6
e_mesh = 0.05
d = 2.1
s = 0.45

es = np.arange(e_min, e_max, e_mesh, dtype=np.complex)

bs  = (es - d) / s
rts = np.sqrt(bs * bs - 4)
z1 = (-bs + rts) / 2
z2 = (-bs - rts) / 2

运行时会产生

z1
Out[118]: 
array([ 0.22222222-0.97499604j,  0.16666667-0.9860133j ,
        0.11111111-0.99380799j,  0.05555556-0.9984556j ,
        0.00000000+1.j        , -0.05555556+0.9984556j ,
       -0.11111111+0.99380799j, -0.16666667+0.9860133j ,
       -0.22222222+0.97499604j, -0.27777778+0.96064536j,
       -0.33333333+0.94280904j, -0.38888889+0.92128466j,
       -0.44444444+0.89580642j, -0.50000000+0.8660254j ,
       -0.55555556+0.83147942j])

z2
Out[119]: 
array([ 0.22222222+0.97499604j,  0.16666667+0.9860133j ,
        0.11111111+0.99380799j,  0.05555556+0.9984556j ,
       -0.00000000-1.j        , -0.05555556-0.9984556j ,
       -0.11111111-0.99380799j, -0.16666667-0.9860133j ,
       -0.22222222-0.97499604j, -0.27777778-0.96064536j,
       -0.33333333-0.94280904j, -0.38888889-0.92128466j,
       -0.44444444-0.89580642j, -0.50000000-0.8660254j ,
       -0.55555556-0.83147942j])