# Large set of nonlinear equations in Sympy

I have a set of 6 nonlinear equations, and using Sympy I find the values of the 6 unknowns. This works perfectly and it directly gives the exact solution, using sympy.solve to be specific. Now I actually have a minimum of 26 equations, and the issue is that once I try to solve for above 6 equations, like 11 equations let us say, then the solutions takes a very long time. I confirmed that my equations are correctly represented so I am not sure how can I use methods in Sympy to check where things are going wrong. I link below a question that I found close to my situation, however I could not benefit from the conclusion of that question:

Speeding up the solution of a large set of nonlinear algebraic equations in sympy

An example I work on for 11 equations is as follows:

from sympy import symbols, I, solve
x1 = symbols('x1', real = True, positive = True)
y1, y2, z1, z2, za1, za2, zb1, zb2, zmj, zjm  = symbols('y1, y2, z1, z2, za1, za2, zb1, zb2, zmj, zjm')
all_symbols = [x, y1, y2, z1, z2, za1, za2, zb1, zb2, zmj, zjm]

eq1 = -x + 0.002*I*(-z1 + y1) + 0.002*I*(-z2 + y2)
eq2 = -0.004*I*x*za1 + 0.002*I*x - 0.5255*y1 - 0.002*I*zmj - 0.002*I*za1
eq3 = -0.004*I*x*za2 + 0.002*I*x + y2*(-0.5255 - 0.5*I) - 0.002*I*zjm - 0.002*I*za2
eq4 = 0.004*I*x*za1 - 0.002*I*x - 0.5255*z1 + 0.002*I*zjm + 0.002*I*za1
eq5 = 0.004*I*x*za2 - 0.002*I*x + z2*(-0.5255 + 0.5*I) + 0.002*I*zmj + 0.002*I*za2
eq6 = 0.002*I*z1 - 0.002*I*y1 - 0.051*za1 + 0.05
eq7 = 0.002*I*z2 - 0.002*I*y2 - 0.051*za2 + 0.05
eq8 = 0.002*I*z1*(1 - 2*za1) - 0.002*I*y1*(1 - 2*za1) - 0.051*zb1
eq9 = 0.002*I*z2*(1 - 2*za2) - 0.002*I*y2*(1 - 2*za2) - 0.051*zb2
eq10 = 0.002*I*z2*(1 - 2*za1) - 0.002*I*y1*(1 - 2*zb2) - 0.051*zmj + 0.5*I*zmj
eq11 = 0.002*I*z1*(1 - 2*za2) - 0.002*I*y2*(1 - 2*zb1) - 0.051*zjm - 0.5*I*zjm

all_equations = [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11]
sol = solve(all_equations, all_symbols, simplify=False)

• As a first approach you can set check=False to speed up solve() but comes at a risk of obtaining invalid solutions. You can also set simplify=False. Depending on your equations it might be worth it to try scipy.optimize.fsolve. Commented Nov 29, 2023 at 11:43
• @mmikkelsen thank you for your comment. I gave the option simplify= False an option, however nothing changed. And for the system I am working with I prefer to have an exact solution and not have to give a first estimate. Commented Nov 29, 2023 at 12:53
• I modified my above text to include part of the code that is in question. Commented Nov 29, 2023 at 13:07
• As you have floating point constants in the equations, you will most likely not get an exact solution, floating-point artifacts will proliferate, terms that should cancel will leave a remainder etc. Multiply the equations by 1000 and change the floating-point numbers to integers. Commented Nov 29, 2023 at 15:19

I tried rewriting your code to use scipy.optimize.fsolve.

import numpy as np
from scipy.optimize import fsolve

def equations(vars):
x, y1, y2, z1, z2, za1, za2, zb1, zb2, zmj, zjm = vars[:11]
I = 1j #To match your notation

eq1 = -x + 0.002*I*(-z1 + y1) + 0.002*I*(-z2 + y2)
eq2 = -0.004*I*x*za1 + 0.002*I*x - 0.5255*y1 - 0.002*I*zmj - 0.002*I*za1
eq3 = -0.004*I*x*za2 + 0.002*I*x + y2*(-0.5255 - 0.5*I) - 0.002*I*zjm - 0.002*I*za2
eq4 = 0.004*I*x*za1 - 0.002*I*x - 0.5255*z1 + 0.002*I*zjm + 0.002*I*za1
eq5 = 0.004*I*x*za2 - 0.002*I*x + z2*(-0.5255 + 0.5*I) + 0.002*I*zmj + 0.002*I*za2
eq6 = 0.002*I*z1 - 0.002*I*y1 - 0.051*za1 + 0.05
eq7 = 0.002*I*z2 - 0.002*I*y2 - 0.051*za2 + 0.05
eq8 = 0.002*I*z1*(1 - 2*za1) - 0.002*I*y1*(1 - 2*za1) - 0.051*zb1
eq9 = 0.002*I*z2*(1 - 2*za2) - 0.002*I*y2*(1 - 2*za2) - 0.051*zb2
eq10 = 0.002*I*z2*(1 - 2*za1) - 0.002*I*y1*(1 - 2*zb2) - 0.051*zmj + 0.5*I*zmj
eq11 = 0.002*I*z1*(1 - 2*za2) - 0.002*I*y2*(1 - 2*zb1) - 0.051*zjm - 0.5*I*zjm

eqs = [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11]
real_parts = [np.real(eq) for eq in eqs]
imag_parts = [np.imag(eq) for eq in eqs]

return real_parts + imag_parts

initial_guess = np.zeros(22)

solution = fsolve(equations, initial_guess)

print(solution)



I had to split the functions into real and imaginary parts since fsolve() only handles real-valued functions. The script above returns

[-1.14007415e-05 -8.51811750e-08 -1.86023150e-03  8.20597016e-08
-1.86029185e-03  9.77386179e-01  9.78812147e-01  3.05561786e-08
-5.27139661e-17 -3.06421661e-05 -1.62352563e-05 -1.39994297e+03
-2.30373227e+02  3.25367804e+02 -1.97975617e+03 -1.35508900e+03
1.62146151e-02  1.46395087e-01 -4.44338695e+02 -4.68748698e+01
9.56867957e+00 -3.68274254e+02]


and when I time it I get 0,44s user 0,29s system 399% cpu 0,183 total. To test the solutions you can write something like this

def test_solution(solution):
x, y1, y2, z1, z2, za1, za2, zb1, zb2, zmj, zjm = solution[:11]
I = 1j

eq1 = -x + 0.002*I*(-z1 + y1) + 0.002*I*(-z2 + y2)
eq2 = -0.004*I*x*za1 + 0.002*I*x - 0.5255*y1 - 0.002*I*zmj - 0.002*I*za1
eq3 = -0.004*I*x*za2 + 0.002*I*x + y2*(-0.5255 - 0.5*I) - 0.002*I*zjm - 0.002*I*za2
eq4 = 0.004*I*x*za1 - 0.002*I*x - 0.5255*z1 + 0.002*I*zjm + 0.002*I*za1
eq5 = 0.004*I*x*za2 - 0.002*I*x + z2*(-0.5255 + 0.5*I) + 0.002*I*zmj + 0.002*I*za2
eq6 = 0.002*I*z1 - 0.002*I*y1 - 0.051*za1 + 0.05
eq7 = 0.002*I*z2 - 0.002*I*y2 - 0.051*za2 + 0.05
eq8 = 0.002*I*z1*(1 - 2*za1) - 0.002*I*y1*(1 - 2*za1) - 0.051*zb1
eq9 = 0.002*I*z2*(1 - 2*za2) - 0.002*I*y2*(1 - 2*za2) - 0.051*zb2
eq10 = 0.002*I*z2*(1 - 2*za1) - 0.002*I*y1*(1 - 2*zb2) - 0.051*zmj + 0.5*I*zmj
eq11 = 0.002*I*z1*(1 - 2*za2) - 0.002*I*y2*(1 - 2*zb1) - 0.051*zjm - 0.5*I*zjm

eqs = [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11]

for i, eq in enumerate(eqs, 1):
print(f"Equation {i}: {eq} (Magnitude: {abs(eq)})")

test_solution(solution)


which will print

Equation 1: (1.1400741530335398e-05-2.1379344165987158e-10j) (Magnitude: 1.140074153233999e-05)
Equation 2: (4.4762707459211876e-08-0.0019546893028441117j) (Magnitude: 0.0019546893033566483)
Equation 3: (0.0009775516553196715-0.0010274542367839554j) (Magnitude: 0.0014181923168257295)
Equation 4: (-4.3122373167240505e-08+0.0019547181166636417j) (Magnitude: 0.0019547181171392955)
Equation 5: (0.0009775833661735019+0.0010273952508865658j) (Magnitude: 0.001418171442161837)
Equation 6: (0.0001533048862500594+3.344817530978207e-10j) (Magnitude: 0.00015330488625042428)
Equation 7: (8.058048974233373e-05-1.206883114378e-10j) (Magnitude: 8.058048974242411e-05)
Equation 8: (-1.558365108814448e-09-3.193539319131462e-10j) (Magnitude: 1.5907510006907593e-09)
Equation 9: (2.6884122732743115e-18+1.1557405909723962e-10j) (Magnitude: 1.1557405909723966e-10)
Equation 10: (1.5627504688856028e-06-1.1768602199419145e-05j) (Magnitude: 1.1871907460731628e-05)
Equation 11: (8.279980708669415e-07+1.1837933761537407e-05j) (Magnitude: 1.18668553773907e-05)


You can properly get better results by tweaking the initial guess and other parameters, here's a link

• thank you for your demonstration with scipy.optimize.fsolve. However I tested the method for the case of 5 equations which I can easily obtain with sympy.solve and I have x = 1.49208588848526e-5, y = -0.00373021472121315*I, z = 0.00373021472121315*I, za = 0.980099591002258, zb = 0.000280921499922358, on the other hand with the fsolve function I have x = 5.01166371e-22 , y = 9.77385972e-01*1j, z = 0.00000000e+00, za = 8.83883388e+01 - 8.83883388e+01*1j, zb = 0.00000000e+00, Since it is very important for me to have the exact solution I can not work with fsolve. Commented Nov 29, 2023 at 14:23
• Have you tried changing the parameters in fsolve? It is very sensitive to initial parameters. You can also use optimize.fmin() to extract a parameter vector. I'm down to approx. 3e-17 for some equations now. Commented Nov 29, 2023 at 15:04
• If you transform a 11-dimensional system for complex variables to real variables, you get a system of 22 variables, with all the difficulties that entails. As the complex root set has a plane for any root of the original system, the "landscape" is extra difficult so finding the real roots could be especially hard. Commented Nov 29, 2023 at 15:30
• @mmikkelsen thanks for the explanation, yes I did try to change the parameters in fsolve, the thing is that first I have complex values in my system and second when I increase the number of equations and change some parameters I can no longer give a close enough estimate so it makes things less efficient, as the exact solution is required. Now for optimize.fmin() I never worked with this before so perhaps I can give it a try, however I am interested in checking if there is a way to make the code work with sympy.solve. Commented Nov 29, 2023 at 17:04
• If I see this right, you have 8 quadratic equations in the system. This makes up to $2^8=256$ roots. I am not aware that sympy implements a Gröbner algorithm, and even then, everything above 50 complex roots is hard. Try a proper CAS, like with python bindings of maxima or directly sage. But still, floating point and symbolic algebra do not work together. Commented Nov 29, 2023 at 18:44