Scipy Fsolve 在具有解的非线性方程组上失败 [英] Scipy Fsolve fails on system of nonlinear equations that has a solution

问题描述

我有一个由 12 个非线性方程组成的系统，如下所示:

I have a system of 12 nonlinear equations as follows:

import math
import numpy as np

# First: Defining some variables that have been calculated before (for the code to work)

n_s1=[8.71557427e-02, 9.96194698e-01, 6.12323400e-17]
n_s2=[5.23359562e-02, 9.98629535e-01, 6.12323400e-17]
P_sens1=[ 2.00000000e+00,  5.01223926e-01, -1.43937571e-17]
P_sens2=[ 2.00000000e+00,  4.20537892e-01, -2.17255041e-17]

# My 12 unknowns, initialized to a known solution: 
ep1_x,ep1_y,ep2_x,ep2_y,er1_x,er1_y,er2_x,er2_y,ep_x,ep_y,en_x,en_y = [0.0287630388891315, 
0.32876303888913166, 0.016591877224378986, 0.3165918772243792, 0.9961946980917455, 0.08715574274765804,
 0.9986295347545739, 0.052335956242943855, 0.0287630388891315, 0.32876303888913166, -0.7071067811865475,
 0.7071067811865476]
# My 12 equations: 
eq = 12*[0]
eq[0] = (ep_x*en_x+ep_y*en_y)/(en_x*n_s1[0]+en_y*n_s1[1])*n_s1[0]-ep1_x
eq[1] = (ep_x*en_x+ep_y*en_y)/(en_x*n_s1[0]+en_y*n_s1[1])*n_s1[1]-ep1_y
eq[2] = n_s1[0]-2*(en_x*n_s1[0]+en_y*n_s1[1])/np.linalg.norm(n_s1)*en_x-er1_x
eq[3] = n_s1[1]-2*(en_x*n_s1[0]+en_y*n_s1[1])/np.linalg.norm(n_s1)*en_y-er1_y
eq[4] = (P_sens1[0]-ep1_x)/math.sqrt(math.pow((P_sens1[0]-ep1_x),2)+math.pow((P_sens1[1]-ep1_y),2))-er1_x
eq[5] = (P_sens1[1]-ep1_y)/math.sqrt(math.pow((P_sens1[0]-ep1_x),2)+math.pow((P_sens1[1]-ep1_y),2))-er1_y
eq[6] = (ep_x*en_x+ep_y*en_y)/(en_x*n_s2[0]+en_y*n_s2[1])*n_s2[0]-ep2_x
eq[7] = (ep_x*en_x+ep_y*en_y)/(en_x*n_s2[0]+en_y*n_s2[1])*n_s2[1]-ep2_y
eq[8] = n_s2[0]-2*(en_x*n_s2[0]+en_y*n_s2[1])/np.linalg.norm(n_s2)*en_x-er2_x
eq[9] = n_s2[1]-2*(en_x*n_s2[0]+en_y*n_s2[1])/np.linalg.norm(n_s2)*en_y-er2_y
eq[10] = (P_sens2[0]-ep2_x)/math.sqrt(math.pow((P_sens2[0]-ep2_x),2)+math.pow((P_sens2[1]-ep2_y),2))-er2_x
eq[11] = (P_sens2[1]-ep2_y)/math.sqrt(math.pow((P_sens2[0]-ep2_x),2)+math.pow((P_sens2[1]-ep2_y),2))-er2_y

我知道系统有一个解决方案，因为当我在这种情况下打印 eq 列表时，它对所有条目的计算结果为零(即最小化).

I know that the system has a solution since when I print the eq list in this case it evaluates to zero for all entries (i.e. being minimized).

print(np.round(eq,9))

>> [-0. -0. -0. -0.  0. -0. -0. -0.  0.  0.  0. -0.]

但是，使用如下实现的scipy fsolve 方法解决系统对我不起作用

However, solving the system using the scipy fsolve method implemented as follows does not work for me

from scipy.optimize import fsolve

def equations(p):
    ep1_x,ep1_y,ep2_x,ep2_y,er1_x,er1_y,er2_x,er2_y,ep_x,ep_y,en_x,en_y = p
    return (eq)

print(fsolve(equations, 12*[0.5]))

我收到一条错误消息阅读

I get an error message reading

>> RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last ten iterations.
  warnings.warn(msg, RuntimeWarning)

我曾尝试更改 fsolve 的参数，但没有任何效果.任何人都可以帮助我，以便求解器可以找到我正在寻找的解决方案吗?干杯.

I have tried changing the parameters of fsolve with no effect. Can anybody help me out such that the solver can find the solution I am looking for? Cheers.

推荐答案

您没有针对 p 点计算方程，因为您在 equation 之外定义了方程> 功能.

You don't evaluate your equations for the point p, since you defined the equations outside the equation function.

将函数更改为

def equations(p):
    eq = np.zeros(p.shape[0])
    ep1_x,ep1_y,ep2_x,ep2_y,er1_x,er1_y,er2_x,er2_y,ep_x,ep_y,en_x,en_y = p
    eq[0] = (ep_x*en_x+ep_y*en_y)/(en_x*n_s1[0]+en_y*n_s1[1])*n_s1[0]-ep1_x
    eq[1] = (ep_x*en_x+ep_y*en_y)/(en_x*n_s1[0]+en_y*n_s1[1])*n_s1[1]-ep1_y
    eq[2] = n_s1[0]-2*(en_x*n_s1[0]+en_y*n_s1[1])/np.linalg.norm(n_s1)*en_x-er1_x
    eq[3] = n_s1[1]-2*(en_x*n_s1[0]+en_y*n_s1[1])/np.linalg.norm(n_s1)*en_y-er1_y
    eq[4] = (P_sens1[0]-ep1_x)/math.sqrt(math.pow((P_sens1[0]-ep1_x),2)+math.pow((P_sens1[1]-ep1_y),2))-er1_x
    eq[5] = (P_sens1[1]-ep1_y)/math.sqrt(math.pow((P_sens1[0]-ep1_x),2)+math.pow((P_sens1[1]-ep1_y),2))-er1_y
    eq[6] = (ep_x*en_x+ep_y*en_y)/(en_x*n_s2[0]+en_y*n_s2[1])*n_s2[0]-ep2_x
    eq[7] = (ep_x*en_x+ep_y*en_y)/(en_x*n_s2[0]+en_y*n_s2[1])*n_s2[1]-ep2_y
    eq[8] = n_s2[0]-2*(en_x*n_s2[0]+en_y*n_s2[1])/np.linalg.norm(n_s2)*en_x-er2_x
    eq[9] = n_s2[1]-2*(en_x*n_s2[0]+en_y*n_s2[1])/np.linalg.norm(n_s2)*en_y-er2_y
    eq[10] = (P_sens2[0]-ep2_x)/math.sqrt(math.pow((P_sens2[0]-ep2_x),2)+math.pow((P_sens2[1]-ep2_y),2))-er2_x
    eq[11] = (P_sens2[1]-ep2_y)/math.sqrt(math.pow((P_sens2[0]-ep2_x),2)+math.pow((P_sens2[1]-ep2_y),2))-er2_y
    return eq

并使用初始点(0.1, ..., 0.1)，即

sol = fsolve(equations, x0=0.1*np.ones(12))

给出具有客观价值的解决方案

gives a solution with objective value

array([-2.01876016e-12,  2.43202680e-11, -4.70620209e-11,  4.23225760e-11,
        4.71587214e-11, -4.30741137e-11,  1.50569834e-12, -2.75213741e-11,
       -2.21203056e-11, -5.79871359e-11,  2.21933583e-11,  5.72545206e-11])

请注意，警告消息仍然存在，只是表明在最近十次迭代中没有取得良好进展.为了获得更好的收敛，您可以将 jacobian 传递给 fsolve.

Note that the warning message remains and just indicates that there was no good progress in the last ten iterations. In order to get better convergence you can pass the jacobian to fsolve.

这篇关于Scipy Fsolve 在具有解的非线性方程组上失败的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！