在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) 在复平面的第四象限中的复杂根，并使用 scipy.optimize 中的 newton() 函数来实现，然后把它们存储到一个一维数组里。我想到的最好方法是通过在第四象限的有限区域内，按照均匀间隔使用 newton() 函数进行暴力搜索，检查找到的根是否是之前找到的重复根，确认这个根确实是一个根，然后把它存入数组中。一旦这个算法完成，我想按照实部从小到大的顺序对数组进行排序。我的问题是：

(i) 我可以创建一个长度不固定的数组吗？这样我就可以在找到值时不断添加进去。

(ii) 我能以某种方式绘制这个函数，以便可视化这些根吗？数学上说它们都在复平面的同一张表面上。

(iii) 有没有更好的方法来找到这些根？我觉得用我的方法会漏掉很多根。

Answer 1

@DrV的回答看起来不错，所以我在这里只想为你问题的第二部分提供另一种思路。一个很有用的方式来可视化复函数的根，就是绘制实部和虚部为零的等高线。也就是说，你可以在一个比较密集的网格上计算z = Dm(...)，然后使用matplotlib的contour函数来绘制z.real为零和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')

这是生成的图。每条红线和蓝线的交点就是一个根。

这个图表明，你不应该期待找到虚部很大的负根。你可以考虑研究贝塞尔函数和汉克尔函数在大参数下的渐进行为，来寻找证明这一点的思路。

Answer 2

一些回答：

(i) 使用列表。数组的大小是固定的，而列表可以随意增加内容。当你往列表里添加一个新的根时，可以检查之前的根是否已经在列表中，比如通过计算 np.amin(np.abs(np.array(a)-b))，其中 a 是已有根的列表，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 有很多选项。如果你只想绘制零点，可以在 countour 的调用中把数字5（轮廓的数量）替换为 [0]（只绘制列出的轮廓）。

当然，你会把我的 (1-C^2) 替换为你自己的函数。唯一需要注意的是，如果函数接收一个复杂数的数组，它必须返回一个形状相同的结果数组，逐点计算。Imshow 需要接收一个标量数组。更多信息，请查看 imshow 的文档。

(iii) 可能会有，但没有通用的方法可以找到任意函数的所有最小值/最大值/零点。（这个函数甚至可能有无限多个根。）你先绘制函数的想法很好。这样你会更容易理解它的行为。