python - 在 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)

我希望使用 scipy.optimize 中的 newton() 找到尽可能多位于复平面第四象限的 Dm(x) 复数根，然后将它们存储到一维数组中。我能想到的最好方法是通过在第四象限的有限部分上以规则的间隔使用 newton() 来暴力破解它，检查根是否是前一个根的重复项，检查根是否确实是一个root，然后将其存储到数组中。一旦该算法完成，我想通过增加实际组件对数组进行排序。我的问题是:

(i) 我可以创建一个长度未定义的数组，以便在找到值时可以继续向其中添加值吗？

(ii) 我可以用可以可视化根的方式绘制函数吗？数学表明它们都在复平面的一张纸上。

(iii) 有没有更好的求根方法？我觉得用我的方法我会错过该领域的很多根源。

Answer 1

一些答案​​:

(i) 使用列表。数组具有固定大小。附加到列表是一个非常便宜的选择。当您将新的根添加到列表中时，请检查前一个根是否不在列表中，例如计算 np.amin(np.abs(np.array(a)-b))哪里a是现 Root过的列表和 b是新的根。如果该值非常小，则表明您已到达现有根。 (多小取决于函数。它不能为 0.0，因为这样您通常会由于浮点和迭代不准确而无法识别相同的根。)

如果您有大量根(数千个)，您可能希望在收到它们后立即对它们进行排序。这使得搜索匹配的根更快。另一方面，很可能超过 90% 的时间都花在迭代根上，您无需担心其他性能问题。然后，您只需编译列表，对其进行排序(列表排序既简单又快速)，然后根据需要转换为数组。

(ii) 是的。下面是两个例子:(对于 countour 的东西，谢谢 Warren 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 :

请注意，如果要缩放图像，通常需要更改限制并重新计算复数值 z以避免遗漏细节。图像当然可以绘制在彼此之上，图像调色板可以更改，countour有一百万个选择。如果您只想绘制零，请将数字 5(轮廓数)替换为 [0] (仅绘制列出的轮廓)在countour中来电。

当然，你可以用你自己的函数替换我的(1-C^2)。唯一需要注意的是，如果函数接收复数数组，它将返回逐点计算的相同形状的结果数组。 Imshow 需要接收一个标量数组。欲了解更多信息，请参阅 imshow文档。

(iii) 可能有，但没有通用方法来查找任意函数的所有最小值/最大值/零点。 (该函数甚至可能有无限多个根。)您首先绘制该函数的想法是一个很好的想法。然后你会更容易理解它的行为。