Scipy Fsolve fails on system of nonlinear equations that has a solution

Question

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

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.]

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)

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.

Answer 1

I solved the problem using python's scipy.optimize.leastsq. In detail the code looks as follows. Also, I can use it for overdetermined systems which I might have in the future:

def f(x):
    ep1_x,ep1_y,ep2_x,ep2_y,er1_x,er1_y,er2_x,er2_y,ep_x,ep_y,en_x,en_y = x
    return np.asarray(((ep_x*en_x+ep_y*en_y)/(en_x*n_s1[0]+en_y*n_s1[1])*n_s1[0]-ep1_x,
                       (ep_x*en_x+ep_y*en_y)/(en_x*n_s1[0]+en_y*n_s1[1])*n_s1[1]-ep1_y,
                       n_s1[0]-2*(en_x*n_s1[0]+en_y*n_s1[1])/np.linalg.norm(n_s1)*en_x-er1_x,
                       n_s1[1]-2*(en_x*n_s1[0]+en_y*n_s1[1])/np.linalg.norm(n_s1)*en_y-er1_y,
                       (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,
                       (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,
                       (ep_x*en_x+ep_y*en_y)/(en_x*n_s2[0]+en_y*n_s2[1])*n_s2[0]-ep2_x,
                       (ep_x*en_x+ep_y*en_y)/(en_x*n_s2[0]+en_y*n_s2[1])*n_s2[1]-ep2_y,
                       n_s2[0]-2*(en_x*n_s2[0]+en_y*n_s2[1])/np.linalg.norm(n_s2)*en_x-er2_x,
                       n_s2[1]-2*(en_x*n_s2[0]+en_y*n_s2[1])/np.linalg.norm(n_s2)*en_y-er2_y,
                       (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,
                       (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))

def system(x,b):
    return (f(x)-b)

b = np.zeros(12)
solution = optimize.leastsq(system, np.ones(12), args=b)[0]

Answer 2

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

Changing the function to

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

and using the initial point (0.1, ..., 0.1), i.e.

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])

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.