存儲複雜功能的根在一個陣列中SciPy的

Question

我有下面的複變函數：

import numpy as np 
import scipy as sp 
from scipy.special import jv, hankel1, jvp, h1vp, h2vp, gamma 

nr = 2 
m = 20 

def Dm(x): 
    return nr * jvp(m,nr*x,1) * hankel1(m,x) - jv(m,nr*x) * h1vp(m,x,1) 

我希望找到DM的許多複雜的根（x）中位於第4象限複雜的平面，因爲我可以從scipy.optimize使用newton（），然後將它們存儲到一維數組中。我可以想到的最好方法是在第四象限的有限部分上以規則間隔使用newton（）來強制它，檢查根是否是前一個根的副本，檢查根是否確實是根，然後將其存儲到數組中。一旦算法完成，我想通過增加真實的組件來對數組進行排序。我的問題是：

（i）我可以創建一個未定義長度的數組，我可以繼續添加值，因爲我發現他們？ （ii）我可以用這樣一種方式繪製函數，以便我可以查看根源嗎？數學說，他們都在複雜的飛機上。 （iii）是否有更好的方法來尋找根源？我覺得我會用我的方法錯過許多領域的根源。

Answer 1

一些答案：

（i）使用一個列表。數組具有固定大小。追加列表是一個非常便宜的選擇。在將新根添加到列表中時，通過（例如）計算np.amin(np.abs(np.array(a)-b))來檢查前一個根是否不在列表中，其中a是現有根的列表，而b是新的根。如果此值非常小，則您已到達現有根目錄。 （如何小取決於函數，它不能爲0.0，因爲通常由於浮點和迭代不準確而無法識別相同的根。）

如果您有大量的根（數千個），那麼您可能希望在收到它們後儘快對其進行分類。這使得更快地搜索匹配的根。另一方面，大概90％以上的時間花在迭代根上，並且不需要擔心其他性能問題。然後，您只需編譯列表，對其進行排序（列表排序非常簡單快捷），並根據需要轉換爲數組。 （ii）是的。下面兩個例子：（對於countour東西，謝謝屬於沃倫Weckesser和他很好的答案！）

import numpy as np 
from scipy.special import jv, hankel1, jvp, h1vp 
import matplotlib.pyplot as plt 


nr = 2 
m = 20 

# create a 2000 x 2000 sample complex plane between -2-2i .. +2+2i 
x = np.linspace(-2, 2, 2000) 
y = np.linspace(-2, 2, 2000) 
X, Y = np.meshgrid(x, y) 
C = X + 1j * Y 
z = 1-C**2 

# draw a contour image of the imaginary part (red) and real part (blue) 
csr = plt.contour(x, y, z.real, 5, colors='b') 
plt.clabel(csr) 
csi = plt.contour(x, y, z.imag, 5, colors='r') 
plt.clabel(csi) 
plt.axis('equal') 
plt.savefig('contours.png') 

# draw an image of the absolute value of the function, black representing zeros 
plt.figure() 
plt.imshow(abs(z), extent=[-2,2,-2,2], cmap=plt.cm.gray) 
plt.axis('equal') 
plt.savefig('absval.png') 

這給countours.png：

和absval.png：

請注意，如果你想放大im年齡時，通常需要更改限制並重新計算複雜值z以避免遺漏細節。這些圖像當然可以相互疊加繪製，圖像調色板可以更改，countour有一百萬個選項。如果您只想繪製零點，請在countour調用中用[0]（僅繪製所列輪廓）替換編號5（輪廓線數量）。

當然，你會用你自己的函數替換我的（1-C^2）。唯一需要注意的是，如果函數接收到一個複數的數組，它將返回一個結果數組，並以逐點形式計算出相同的形狀。 Imshow需要接收一組標量。有關更多信息，請參閱imshow文檔。 （iii）可能有，但沒有找到任意函數的所有最小/最大/零的一般方法。 （該函數甚至可能有無數的根。）您首先繪製函數的想法是一個好主意。然後你就會更容易理解它的行爲。

Answer 2

@ DrV的答案看起來不錯，所以在這裏我只會提出另一種方法來處理你的問題的（ii） 。可視化複雜函數根的有用方法是繪製實部和虛部的0輪廓。也就是說，計算 z = Dm(...)在一個合理的密網格，然後用matplotlib的contour功能 畫出輪廓，其中z.real爲0，其中z.imag爲零。函數的根是這些輪廓相交的點。

例如，下面是一些生成輪廓圖的代碼。

import numpy as np 
from scipy.special import jv, hankel1, jvp, h1vp 
import matplotlib.pyplot as plt 


nr = 2 
m = 20 

def Dm(x): 
    return nr * jvp(m,nr*x,1) * hankel1(m,x) - jv(m,nr*x) * h1vp(m,x,1) 


x = np.linspace(0, 40, 2500) 
y = np.linspace(-15, 5, 2000) 
X, Y = np.meshgrid(x, y) 
z = Dm(X + 1j*Y) 

plt.contour(x, y, z.real, [0], colors='b', lw=0.25) 
plt.contour(x, y, z.imag, [0], colors='r', lw=0.25) 
plt.savefig('contours.png') 

這裏是情節。紅線和藍線的每個交點都是根。

情節表明，你不應該期望找到與大 負虛部的根。你或許可以看看貝塞爾函數和漢克爾函數的漸近行爲 ，其中有很多想法可以證明這一點。