在复平面和corr中找到A(x)和B(y)的交点.x和y

Question

假设我有以下问题:

我有一个复杂的功能A(x)和复杂的功能B(y).我知道这些函数在复杂的平面上交叉.我想用数字(和/或图形)找出相应的x和y这个交叉点.这样做最聪明的方法是什么？

这是我的出发点:

import matplotlib.pyplot as plt
import numpy as np
from numpy import sqrt, pi

x = np.linspace(1, 10, 10000)
y = np.linspace(1, 60, 10000)

def A_(x):
    return -1/( 8/(pi*x)*sqrt(1-(1/x)**2) - 1j*(8/(pi*x**2)) )
A = np.vectorize(A_)

def B_(y):
    return 3/(1j*y*(1+1j*y))
B = np.vectorize(B_)


real_A = np.real(A(x))
imag_A = np.imag(A(x))

real_B = np.real(B(y))
imag_B = np.imag(B(y))

plt.plot(real_A, imag_A, color='blue')
plt.plot(real_B, imag_B, color='red')
plt.show()

我不一定要画它.我只需要x_intersection和y_intersection(有一些错误取决于x和y).

非常感谢提前!

编辑: 我应该使用不同的变量名称.为了澄清我需要: x和y是numpy的阵列和我需要每个阵列的交点的索引加上相应的x和y值(这又是不相交点本身,而是阵列的一些值x和y).

Answer 1

在这里,我找到两条曲线之间的距离的最小值.此外,我稍微清理了你的代码(例如,vectorize没有做任何有用的事情).

import matplotlib.pyplot as plt
import numpy as np
from numpy import sqrt, pi
from scipy import optimize

def A(x):
    return -1/( 8/(pi*x)*sqrt(1-(1/x)**2) - 1j*(8/(pi*x**2)) )

def B(y):
    return 3/(1j*y*(1+1j*y))

#  The next three lines find the intersection
def dist(x):
    return abs(A(x[0])-B(x[1]))    
sln = optimize.minimize(dist, [1, 1])

#  plotting everything....
a0, b0 = A(sln.x[0]), B(sln.x[1])
x = np.linspace(1, 10, 10000)
y = np.linspace(1, 60, 10000)    
a, b = A(x), B(y)

plt.plot(a.real, a.imag, color='blue')
plt.plot(b.real, b.imag, color='red')
plt.plot(a0.real, a0.imag, "ob")
plt.plot(b0.real, b0.imag, "xr")
plt.show()

交点处的特定x和y值是sln.x[0]和sln.x[1],因为A(sln.x[0])=B(sln.x[1]).如果你需要的指数,这也是你在你的编辑提,我会使用,例如numpy.searchsorted(x, sln.x[0]),寻找到从拟合值将插入到您x和y阵列.

我认为这个问题有点棘手的问题是交叉点所在的图形空间(即复平面)不显示输入空间,但必须优化输入空间.它可用于可视化解决方案,然后绘制输入空间上的曲线之间的距离.这可以这样做:

data = dist((X, Y))
fig, ax = plt.subplots()
im = ax.imshow(data, cmap=plt.cm.afmhot, interpolation='none', 
    extent=[min(x), max(x), min(y), max(y)], origin="lower")
cbar = fig.colorbar(im)
plt.plot(sln.x[0], sln.x[1], "xw")
plt.title("abs(A(x)-B(y))")

从这一点来看,似乎更清楚如何optimize.minimum工作 - 它只是在斜坡上滚动以找到最小距离,在这种情况下为零.但是,仍然没有明显的单一可视化可以用来查看整个问题.

对于其他交叉路口,人们必须多挖一点.也就是说,@ emma在评论中询问了其他根源,并且在那里我提到没有一般可靠的方法来找到任意方程式的所有根,但这里是我如何寻找其他根.在这里,我不会列出完整的程序,只是列出了我的变化和情节.

首先,显而易见的是,对于我的第一个图中显示的域,只有一个交叉点,并且左边区域没有交叉点.唯一可能有另一个交叉点的地方是右边的,但是为了那个我需要允许sqrt在def中B获得否定参数而不抛出异常.一个简单的方法是添加0j这样的参数sqrt,像这样sqrt(1+0j-(1/x)**2).然后与交叉点的情节成为

我在更广泛的范围(x=np.linspace(-10, 10, 10000)和y=np.linspace(-400, 400, 10000))上绘制了这个,上面是唯一有趣的地方的缩放.这显示了上面找到的交点,加上看起来两条曲线可能接触的点(红色曲线B到达几乎与蓝色曲线A向上相交的点),这是新的有趣的东西,我的事情我会找.

有点玩限制等,表明B渐渐地达到了一个点,并且B显而易见的是它会变得0 + 0j很大+/- y,所以这就是所有可以说的B.

A从上面的图中很难理解,所以我将独立地看一下实部和虚部:

所以它不是一个看起来很疯狂的函数,而且Re = const和Im = const之间的跳跃只是sqrt(1-x ^-2)的本质,它纯粹是复杂的abs(x)<1,纯粹的真实abs(x)>1.

现在非常清楚,曲线相等的另一个时间是y= +/-inf和x=0.并且,快速查看方程式显示A(0)=0+0j和B(+/- inf)=0+0j,所以这是另一个交叉点(虽然它发生在B(+/- inf),但是它是否会被称为交叉点是不明确的).

这就是它.另一点要提的是,如果这些没有这样一个简单的解析解,就像是不清楚什么B是在inf等,一个也图/最小化等,通过观察B(1/y),然后从那里走,使用与上面相同的工具来处理无穷大.所以使用:

def dist2(x):
    return abs(A(x[0])-B(1./x[1]))

右边的min是最初找到的那个,而零,现在是x=-0和1./y=0另一个(再次,这里没有足够的兴趣应用优化器,但在其他方程中可能很有趣).

当然,也可以通过查找上图中的最小数据来估计这一点,如下所示:

X, Y = np.meshgrid(x, y)
data = dist((X, Y))
r = np.unravel_index(data.argmin(), data.shape)
print x[r[1]], y[r[0]]
#  2.06306306306 1.8008008008    # min approach gave 2.05973231  1.80069353

但这只是近似(对于分辨率data)并且涉及更多的计算(1M与几百相比).我只发布这个,因为我认为这可能是OP最初的想法.