如何求解具有三个未知数和数百个解的 Python 的非线性方程？
How to solve nonlinear equation with Python with three unknowns and hundreds of solutions?
我正在尝试使用 python 在以下类型的非线性方程中求出三个未知数 (x,y,z) 的值：
g(x) * h(y) * k(z) = F
其中 F 是具有数百个值的向量。
我成功地使用了 scipy.optimize.minimize，其中 F 只有 3 个值，但是当 F 的大小大于 3 时失败了。
如何使用 F 中的所有值找到 (x,y,z)？到目前为止，唯一可用的方法是使用网格搜索，但这效率很低（请参阅# Find unknown x[0], x[1], x[2]）。 python 中是否有一个函数可以用来查找 x[0]、x[1]、x[2] 而不是使用网格搜索？
这里是代码：
import numpy as np
# Inputs:
thetas = np.array([25.4,65,37,54.9,26,21.3,24.1,35.7,46.1,61.1,57.2,41.9,20.5,24,55.6,56.9,42.2,39.9,30.8,59,28.8])
thetav = np.array([28.7,5.4,22.6,14.4,23.5,25,12.8,31.2,15.3,9,7.4,24.4,29.7,15.3,15.5,26.8,8.8,16.6,25.1,18.5,12])
azs = np.array([130.3,158,150.2,164.8,152.4,143.5,144.2,151.8,167.4,169.7,162.2,161.4,138.2,147.8,172.9,168.6,158.3,159.8,151.7,160.8,144.5])
azv = np.array([55.9,312.8,38.6,160.4,324.2,314.5,236.3,86.1,313.3,2.1,247.6,260.4,118.9,199.9,277,103.1,150.5,339.2,35.6,14.7,24.9])
F = np.array([0.61745,0.43462,0.60387,0.56595,0.48926,0.55615,0.54351,0.64069,0.54228,0.51716,0.39157,0.51831,0.7053,0.62769,0.21159,0.29964,0.52126,0.53656,0.575,0.40306,0.60471])
relphi = np.abs(azs-azv)
thetas = np.deg2rad(thetas)
thetav = np.deg2rad(thetav)
relphi = np.deg2rad(relphi)
# Compute the trigonometric functions:
coss = np.abs (np.cos(thetas))
cosv = np.cos(thetav)
sins = np.sqrt(1.0 - coss * coss)
sinv = np.sqrt(1.0 - cosv * cosv)
cosp = -np.cos(relphi)
tans = sins / coss
tanv = sinv / cosv
csmllg = coss * cosv + sins * sinv * cosp
bigg = np.sqrt(tans * tans + tanv * tanv - 2.0 * tans * tanv * cosp)
# Function to solve
def fun(x):
return x[0] * ((coss * cosv) ** (x[1] - 1.0)) * ((coss + cosv) ** (x[1] - 1.0)) * (1.0 - x[2] * x[2]) / ((1.0 + x[2] * x[2] + 2.0 * x[2] * csmllg) ** (1.5) + 1e-12) * (1.0 + ((1 - x[0]) / (1.0 + bigg))) - F
# Find unknown x[0], x[1], x[2]
n_bins=51
rho0_min=0.0
rho0_max=2.0
rho0_index=np.linspace(rho0_min, rho0_max, n_bins,retstep=True)
k_min=0.0
k_max=2.0
k_index=np.linspace(k_min, k_max, n_bins,retstep=True)
bigtet_min=-1.0
bigtet_max=1.0
bigtet_index=np.linspace(bigtet_min, bigtet_max, n_bins,retstep=True)
results=np.zeros((4,n_bins**3))
minima=np.ones(4)
RMSE_th = 0.001
index_while=0
current_RMSE=1.0
while current_RMSE > RMSE_th:
index_results=0
for rho0 in rho0_index[0]:
for k in k_index[0]:
for bigtet in bigtet_index[0]:
results[:,index_results] = [rho0, k, bigtet, np.sqrt(np.sum((surf-func([rho0,k,bigtet]))**2) / surf.size)]
index_results=index_results+1
minima = results[:,np.argmin(results[3,:])]
if (index_while > 10) or ((current_RMSE-minima[3]) < RMSE_th/100.0):
break
else:
current_RMSE=minima[3]
index_while=index_while+1
rho0_min=minima[0]-2*rho0_index[1]
rho0_max=minima[0]+2*rho0_index[1]
rho0_index=np.linspace(rho0_min, rho0_max, 11,retstep=True)
k_min=minima[1]-2*k_index[1]
k_max=minima[1]+2*k_index[1]
k_index=np.linspace(k_min, k_max, 11,retstep=True)
bigtet_min=minima[2]-2*bigtet_index[1]
bigtet_max=minima[2]+2*bigtet_index[1]
bigtet_index=np.linspace(bigtet_min, bigtet_max, 11,retstep=True)
rho0= minima[0]
k= minima[1]
bigtet= minima[2]
print (rho0,k,bigtet,minima[3])
return (rho0,k,bigtet)
啊哈，我明白了。您正在使用 scipy.optimize.minimize，它获取函数 fun(x) 作为要最小化的函数。但是根据 documentation 最小化执行 "Minimization of scalar function of one or more variables."。这里的线索是 标量。您需要 return 一个标量。一个经典的方法是最小化方差之和：
deviations = x[0] * ((coss * cosv) ** (x[1] - 1.0)) * ((coss + cosv) ** (x[1] - 1.0)) * (1.0 - x[2] * x[2]) / ((1.0 + x[2] * x[2] + 2.0 * x[2] * csmllg) ** (1.5) + 1e-12) * (1.0 + ((1 - x[0]) / (1.0 + bigg))) - F
我强烈建议您开始将该表达式分解成多个部分，以便更容易理解。那么，函数的 return 可以是：
return sum( deviations ** 2)
正如@Jblasco 所建议的，您可以最小化平方和。 scipy.leastsq()就是为这类问题设计的。对于您的示例，代码为：
import scipy.optimize as sopt
xx0 = np.array([0., 0., 0.]) # starting point
rslt = sopt.leastsq(fun, xx0, full_output=True)
print("The solution is {}".format(rslt[0]))
查看 rslts 的其他条目以获取有关解决方案质量的信息。请记住，数字可能会给您带来麻烦，尤其是在有指数和数百个变量时。如果您遇到问题，请查看 Scipy 中的其他 optimizers。还提供明确的雅可比矩阵（作为 leastsq() 的参数）可以提供帮助。
我正在尝试使用 python 在以下类型的非线性方程中求出三个未知数 (x,y,z) 的值：
g(x) * h(y) * k(z) = F
其中 F 是具有数百个值的向量。
我成功地使用了 scipy.optimize.minimize，其中 F 只有 3 个值，但是当 F 的大小大于 3 时失败了。
如何使用 F 中的所有值找到 (x,y,z)？到目前为止，唯一可用的方法是使用网格搜索，但这效率很低（请参阅# Find unknown x[0], x[1], x[2]）。 python 中是否有一个函数可以用来查找 x[0]、x[1]、x[2] 而不是使用网格搜索？
这里是代码：
import numpy as np
# Inputs:
thetas = np.array([25.4,65,37,54.9,26,21.3,24.1,35.7,46.1,61.1,57.2,41.9,20.5,24,55.6,56.9,42.2,39.9,30.8,59,28.8])
thetav = np.array([28.7,5.4,22.6,14.4,23.5,25,12.8,31.2,15.3,9,7.4,24.4,29.7,15.3,15.5,26.8,8.8,16.6,25.1,18.5,12])
azs = np.array([130.3,158,150.2,164.8,152.4,143.5,144.2,151.8,167.4,169.7,162.2,161.4,138.2,147.8,172.9,168.6,158.3,159.8,151.7,160.8,144.5])
azv = np.array([55.9,312.8,38.6,160.4,324.2,314.5,236.3,86.1,313.3,2.1,247.6,260.4,118.9,199.9,277,103.1,150.5,339.2,35.6,14.7,24.9])
F = np.array([0.61745,0.43462,0.60387,0.56595,0.48926,0.55615,0.54351,0.64069,0.54228,0.51716,0.39157,0.51831,0.7053,0.62769,0.21159,0.29964,0.52126,0.53656,0.575,0.40306,0.60471])
relphi = np.abs(azs-azv)
thetas = np.deg2rad(thetas)
thetav = np.deg2rad(thetav)
relphi = np.deg2rad(relphi)
# Compute the trigonometric functions:
coss = np.abs (np.cos(thetas))
cosv = np.cos(thetav)
sins = np.sqrt(1.0 - coss * coss)
sinv = np.sqrt(1.0 - cosv * cosv)
cosp = -np.cos(relphi)
tans = sins / coss
tanv = sinv / cosv
csmllg = coss * cosv + sins * sinv * cosp
bigg = np.sqrt(tans * tans + tanv * tanv - 2.0 * tans * tanv * cosp)
# Function to solve
def fun(x):
return x[0] * ((coss * cosv) ** (x[1] - 1.0)) * ((coss + cosv) ** (x[1] - 1.0)) * (1.0 - x[2] * x[2]) / ((1.0 + x[2] * x[2] + 2.0 * x[2] * csmllg) ** (1.5) + 1e-12) * (1.0 + ((1 - x[0]) / (1.0 + bigg))) - F
# Find unknown x[0], x[1], x[2]
n_bins=51
rho0_min=0.0
rho0_max=2.0
rho0_index=np.linspace(rho0_min, rho0_max, n_bins,retstep=True)
k_min=0.0
k_max=2.0
k_index=np.linspace(k_min, k_max, n_bins,retstep=True)
bigtet_min=-1.0
bigtet_max=1.0
bigtet_index=np.linspace(bigtet_min, bigtet_max, n_bins,retstep=True)
results=np.zeros((4,n_bins**3))
minima=np.ones(4)
RMSE_th = 0.001
index_while=0
current_RMSE=1.0
while current_RMSE > RMSE_th:
index_results=0
for rho0 in rho0_index[0]:
for k in k_index[0]:
for bigtet in bigtet_index[0]:
results[:,index_results] = [rho0, k, bigtet, np.sqrt(np.sum((surf-func([rho0,k,bigtet]))**2) / surf.size)]
index_results=index_results+1
minima = results[:,np.argmin(results[3,:])]
if (index_while > 10) or ((current_RMSE-minima[3]) < RMSE_th/100.0):
break
else:
current_RMSE=minima[3]
index_while=index_while+1
rho0_min=minima[0]-2*rho0_index[1]
rho0_max=minima[0]+2*rho0_index[1]
rho0_index=np.linspace(rho0_min, rho0_max, 11,retstep=True)
k_min=minima[1]-2*k_index[1]
k_max=minima[1]+2*k_index[1]
k_index=np.linspace(k_min, k_max, 11,retstep=True)
bigtet_min=minima[2]-2*bigtet_index[1]
bigtet_max=minima[2]+2*bigtet_index[1]
bigtet_index=np.linspace(bigtet_min, bigtet_max, 11,retstep=True)
rho0= minima[0]
k= minima[1]
bigtet= minima[2]
print (rho0,k,bigtet,minima[3])
return (rho0,k,bigtet)
啊哈，我明白了。您正在使用 scipy.optimize.minimize，它获取函数 fun(x) 作为要最小化的函数。但是根据 documentation 最小化执行 "Minimization of scalar function of one or more variables."。这里的线索是 标量。您需要 return 一个标量。一个经典的方法是最小化方差之和：
deviations = x[0] * ((coss * cosv) ** (x[1] - 1.0)) * ((coss + cosv) ** (x[1] - 1.0)) * (1.0 - x[2] * x[2]) / ((1.0 + x[2] * x[2] + 2.0 * x[2] * csmllg) ** (1.5) + 1e-12) * (1.0 + ((1 - x[0]) / (1.0 + bigg))) - F
我强烈建议您开始将该表达式分解成多个部分，以便更容易理解。那么，函数的 return 可以是：
return sum( deviations ** 2)
正如@Jblasco 所建议的，您可以最小化平方和。 scipy.leastsq()就是为这类问题设计的。对于您的示例，代码为：
import scipy.optimize as sopt
xx0 = np.array([0., 0., 0.]) # starting point
rslt = sopt.leastsq(fun, xx0, full_output=True)
print("The solution is {}".format(rslt[0]))
查看 rslts 的其他条目以获取有关解决方案质量的信息。请记住，数字可能会给您带来麻烦，尤其是在有指数和数百个变量时。如果您遇到问题，请查看 Scipy 中的其他 optimizers。还提供明确的雅可比矩阵（作为 leastsq() 的参数）可以提供帮助。