尝试找到满足 f(x,y) = 0 和 g(x,y) = 0 的所有有序对 (x,y)；即，我试图找到多变量函数的根
Trying to find all ordered pairs (x,y) such that f(x,y) = 0 and g(x,y) = 0; i.e., I'm trying to find the roots of a multi-variable function
注意：我最初有 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)
这个 returns f 和 g 值数组。
我想找到 f(X,Y) = 0 和 g(X,Y) = 0 的所有 (X,Y) 有序对。我正在寻找不同的 scipy包，我找不到任何似乎适用于这样的多变量函数的东西。另外，我现在的答案是以数组形式返回的，所以我可以使用 np.where() 之类的东西吗？问题在于，因为我存储的是精确值，所以我不一定会看到 f(x,y) 或 g(x,y) 明确等于零。我的最终目标是绘制这些点。此外，到目前为止我所做的对 Xs 和 Ys 作为这些范围内的 linspaces 是否有意义？
谢谢
更新：我回去并使用我在一个有点类似的问题上找到的指南编写了一个小脚本，参见 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) 方面的一些混淆。
如 has already noted in a comment, you probably need scipy.optimize to do the bulk of the work for you. Specifically, either scipy.optimize.fsolve or 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])
wich 将尝试从 [0.5,0.5] 开始寻找零。请注意，向量值函数的根恰好位于其两个分量均为零的位置，这就是您所追求的。
return 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会告诉你函数的值（接近0如果 sol.status==1).
现在，正如您也注意到的那样，每次调用最多只会告诉您一个根。要找到多个根，您无法避免搜索它们。为此，您可以遍历您选择的 X,Y 网格，从那里开始 root/fsolve，然后检查搜索是否成功。如果是：存储 post-processing.
的值
不幸的是，找到非线性多元函数的零点绝非易事，因此您迟早要亲自动手。
更新
你遇到麻烦了。考虑：
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 参数，以确保它们计算相同的东西。
注意：我最初有 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)
这个 returns f 和 g 值数组。
我想找到 f(X,Y) = 0 和 g(X,Y) = 0 的所有 (X,Y) 有序对。我正在寻找不同的 scipy包，我找不到任何似乎适用于这样的多变量函数的东西。另外，我现在的答案是以数组形式返回的，所以我可以使用 np.where() 之类的东西吗？问题在于，因为我存储的是精确值，所以我不一定会看到 f(x,y) 或 g(x,y) 明确等于零。我的最终目标是绘制这些点。此外，到目前为止我所做的对 Xs 和 Ys 作为这些范围内的 linspaces 是否有意义？
谢谢
更新：我回去并使用我在一个有点类似的问题上找到的指南编写了一个小脚本，参见 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) 方面的一些混淆。
如scipy.optimize to do the bulk of the work for you. Specifically, either scipy.optimize.fsolve or 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])
wich 将尝试从 [0.5,0.5] 开始寻找零。请注意，向量值函数的根恰好位于其两个分量均为零的位置，这就是您所追求的。
return 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会告诉你函数的值（接近0如果 sol.status==1).
现在，正如您也注意到的那样，每次调用最多只会告诉您一个根。要找到多个根，您无法避免搜索它们。为此，您可以遍历您选择的 X,Y 网格，从那里开始 root/fsolve，然后检查搜索是否成功。如果是：存储 post-processing.
不幸的是，找到非线性多元函数的零点绝非易事，因此您迟早要亲自动手。
更新
你遇到麻烦了。考虑：
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 参数，以确保它们计算相同的东西。