python - 尝试找到所有满足 f(x,y) = 0 且 g(x,y) = 0 的有序对 (x,y)；即，我试图找到多变量函数的根

Question

注意:我最初有 chi = 1，但我已将其更改为 chi = 0(这是更简单的情况)。

我的方程 f(x,y) 和 g(x,y) 来自以下代码:

import numpy as np
from pylab import *

def stress(X,Y):
    chi = 0
    F = 1
    a = 1
    c = (1.0*a)/(np.sqrt(np.power(X,2)+np.power(Y,2))*1.0)
    A = 0.5*(1 - c**2. + (1 - 4*c**2 + 3*c**4.)*np.cos(2*np.arctan(Y/X)))
    B = 0.5*(1 - c**2. - (1 - 4*c**2 + 3*c**4.)*np.cos(2*np.arctan(Y/X)))
    C = 0.5*(1 + c**2. - (1 + 3*c**4.)*np.cos(2*np.arctan(Y/X)))
    D = 0.5*(1 + c**2. + (1 + 3*c**4.)*np.cos(2*np.arctan(Y/X)))
    E = 0.5*((1 + 2*c**2. - 3*c**4.)*np.sin(2*np.arctan(Y/X)))

    f = 1.0*F*c**2 + (A-1.0*chi*B) # Radial stress
    g = -1.*F*c**2 + (C - 1.0*chi*D) # Tangential stress
    return f,g

def f(X,Y):
    return stress(X,Y)[0]
def g(X,Y):
    return stress(X,Y)[1]

def find_points(X_max,Y_max,X_steps,Y_steps):
    Xs = np.linspace(-X_max,X_max,X_steps)
    Ys = np.linspace(-Y_max,Y_max,Y_steps)

    radials = f(Xs,Ys)
    tangentials = g(Xs,Ys)

    return radials, tangentials

find_points(10,10,100,100)

这将返回 f 和 g 的值数组。

我想找到所有 f(X,Y) = 0 和 g(X,Y) = 0 的 (X,Y) 有序对。我正在查看不同的 scipy 包，但我找不到找到任何似乎适用于这样的多变量函数的东西。另外，我现在的答案正在数组中返回，所以我可以使用像 np.where() 这样的东西吗？问题在于，因为我存储的是精确值，所以我不一定会看到 f(x,y) 或 g(x,y) 明确等于零。我的最终目标是绘制这些点。另外，我到目前为止所做的将 X 和 Y 作为这些范围内的 linspace 有意义吗？

谢谢

更新:我回去并使用我在类似问题上找到的指南编写了一个小脚本，请参阅 This link 。我使用了 scipy.optimize。

from scipy import optimize

def equations(p):
    X, Y = p
    a = 1
    F = 1
    chi = 0
    c = (1.0*a)/(np.sqrt(np.power(X,2)+np.power(Y,2))*1.0)
    A = 0.5*(1 - c**2. + (1 - 4*c**2 + 3*c**4.)*np.cos(2*np.arctan(Y/X)))
    B = 0.5*(1 - c**2. - (1 - 4*c**2 + 3*c**4.)*np.cos(2*np.arctan(Y/X)))
    C = 0.5*(1 + c**2. - (1 + 3*c**4.)*np.cos(2*np.arctan(Y/X)))
    D = 0.5*(1 + c**2. + (1 + 3*c**4.)*np.cos(2*np.arctan(Y/X)))

    f = 1.0*F*c**2 + (A-1.0*chi*B) # Radial stress
    g = -1.*F*c**2 + (C - 1.0*chi*D) # Tangential stress
    return (f,g)

X, Y =  optimize.fsolve(equations, (1, 1))

print equations((X, Y))

这需要我输入不同的初始猜测以获得不同的 (X,Y) 根。如果我能以某种方式解决所有的问题，那就太棒了。另外，我得到的答案似乎有点不对劲。再次感谢。

注意:在将它们转换为笛卡尔坐标之前，原始方程如下:

def stress(R,theta):
    chi = 0
    F = 1
    a = 1
    c = (1.0*a)/(R*1.0)
    A = 0.5*(1 - c**2. + (1 - 4*c**2 + 3*c**4.)*np.cos(2*theta))
    B = 0.5*(1 - c**2. - (1 - 4*c**2 + 3*c**4.)*np.cos(2*theta))
    C = 0.5*(1 + c**2. - (1 + 3*c**4.)*np.cos(2*theta))
    D = 0.5*(1 + c**2. + (1 + 3*c**4.)*np.cos(2*theta))
    E = 0.5*((1 + 2*c**2. - 3*c**4.)*np.sin(2*theta))

    f = 1.0*F*c**2. + (A-1.0*chi*B) # Radial stress
    g = -1.0*F*c**2. + (C-1.0*chi*D) # Tangential stress
    return f,g

也许这将有助于解决方程的 arctan(Y/X) 方面的一些困惑。

Answer 1

如@Azad已经在评论中指出，您可能需要 scipy.optimize为您完成大部分工作。具体来说，scipy.optimize.fsolve或scipy.optimize.root 。由于后者似乎更普遍，我将证明这一点。由于它可以使用多种方法，请查看帮助。

这两个函数都能够找到从 R^n 到 R^m 映射的函数的根，即多元向量值函数。如果您考虑您的stress 函数，这正是您所拥有的:它从 R^2 映射到 R^2。为了清楚起见，您甚至可以将其定义为

def stress2(Rvec):
    X,Y=Rvec
    chi = 1
    F = 1
    a = 1
    c = (1.0*a)/(np.sqrt(np.power(X,2)+np.power(Y,2))*1.0)
    A = 0.5*(1 - c**2. + (1 - 4*c**2 + 3*c**4.)*np.cos(2*np.arctan(Y/X)))
    B = 0.5*(1 - c**2. - (1 - 4*c**2 + 3*c**4.)*np.cos(2*np.arctan(Y/X)))
    C = 0.5*(1 + c**2. - (1 + 3*c**4.)*np.cos(2*np.arctan(Y/X)))
    D = 0.5*(1 + c**2. + (1 + 3*c**4.)*np.cos(2*np.arctan(Y/X)))
    E = 0.5*((1 + 2*c**2. - 3*c**4.)*np.sin(2*np.arctan(Y/X)))
    f = 1.0*F*c**2 + (A-1.0*chi*B) # Radial stress
    g = -1.*F*c**2 + (C - 1.0*chi*D) # Tangential stress
    return f,g

现在，通过这个定义，您可以简单地调用

import scipy.optimize as opt
sol=opt.root(stress2,[0.5,0.5])

它将尝试从[0.5,0.5]开始找到零。请注意，向量值函数的根恰好位于其两个分量均为零的位置，这就是您所追求的。

返回OptimizeResult看起来像这样:

In [224]: sol
Out[224]: 
  status: 1
 success: True
     qtf: array([  2.94481987e-09,   4.76366933e-25])
    nfev: 47
       r: array([ -7.62669534e-06,   7.62669532e-06,   2.16965211e-21])
     fun: array([  2.25125258e-10,  -2.25125258e-10])
       x: array([ 167337.87789902,  167337.87786433])
 message: 'The solution converged.'
    fjac: array([[-0.70710678,  0.70710678],
       [ 0.70710678,  0.70710678]])

它有很多信息。首先，sol.status会告诉你它是否收敛成功。这是最重要的输出:您的根以及非常敏感找到它的可能性取决于您的起点。如果您在示例中尝试从 X=0 或 Y=0 开始，您会发现很难找到根。

如果您确实有根，sol.x会告诉您坐标，sol.fun会告诉您你的函数(如果 sol.status==1 则接近 0)。

现在，正如您还注意到的那样，每个调用最多只会告诉您一个根。要找到多个根，您就无法避免搜索它们。您可以通过检查您选择的 X,Y 网格，从那里启动 root/fsolve 并检查搜索是否成功来完成此操作。如果是:存储值以进行后处理。

不幸的是，找到非线性多元函数的零点远非易事，因此您迟早必须亲自动手。

更新

你遇到麻烦了。考虑:

v=np.linspace(-10,10,100)
X,Y=np.meshgrid(v,v)

fig = plt.figure()
hc=plt.contourf(X,Y,stress2([X,Y])[0].clip(-1,1),levels=np.linspace(-1,1,20))
plt.contour(X,Y,stress2([X,Y])[0].clip(-1,1),levels=[0],color=(1,0,0))
plt.colorbar(hc)

其他函数也是如此。函数的 x 和 y 组件分别如下所示:

它们都沿着一些类似双曲线的曲线有零点。 这似乎是相同的。该图强烈表明，存在一条线点，您的函数为零:两个分量。对于数值求根算法来说，这是最糟糕的，因为没有明确的(孤立的)零。

我建议在纸上检查您的函数是否存在 X==Y 的情况，您可能确实会看到您的函数在那里消失，至少是渐近地消失。

更新2

您添加了函数的原始极坐标形式。虽然我看不出你哪里出了问题(除了使用 np.arctan 而不是 np.arctan2 ，但这似乎并不能解决问题)，我尝试绘制你的极坐标函数:

def stress_polar(Rvec):
    R,theta=Rvec
    chi = 0
    F = 1
    a = 1
    c = (1.0*a)/(R*1.0)                                   
    A = 0.5*(1 - c**2. + (1 - 4*c**2 + 3*c**4.)*np.cos(2*theta))
    B = 0.5*(1 - c**2. - (1 - 4*c**2 + 3*c**4.)*np.cos(2*theta))
    C = 0.5*(1 + c**2. - (1 + 3*c**4.)*np.cos(2*theta))
    D = 0.5*(1 + c**2. + (1 + 3*c**4.)*np.cos(2*theta))
    E = 0.5*((1 + 2*c**2. - 3*c**4.)*np.sin(2*theta))
    f = 1.0*F*c**2. + (A-1.0*chi*B)
    g = -1.0*F*c**2. + (C-1.0*chi*D)
    return f,g

v1=np.linspace(0.01,10,100)
v2=np.linspace(-np.pi,np.pi,100)
R,theta=np.meshgrid(v1,v2)

fig = plt.figure()
ax=plt.subplot(111, polar=True)
hc=plt.contourf(theta,R,stress_polar([R,theta])[0].clip(-1,1),levels=np.linspace(-1,1,20))
plt.contour(theta,R,stress_polar([R,theta])[0].clip(-1,1),levels=[0],color=(1,0,0))
plt.colorbar(hc)

对于切向分量也是如此。请注意，极坐标图需要首先获取 theta，然后获取 R。结果:

这显示了一幅完全不同的图片，对径向分量的零点有有限的支持。现在，我以前从未在 matplotlib 中使用过极坐标图，因此同样有可能我在绘图过程中搞砸了一些东西。但可能值得查看由极坐标函数和笛卡尔函数计算的 A,B,C,D 参数，以确保它们计算相同的东西。