2
$\begingroup$

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)
$\endgroup$
4
  • 1
    $\begingroup$ 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. $\endgroup$
    – mmikkelsen
    Commented Nov 29, 2023 at 11:43
  • $\begingroup$ @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. $\endgroup$
    – je2703
    Commented Nov 29, 2023 at 12:53
  • $\begingroup$ I modified my above text to include part of the code that is in question. $\endgroup$
    – je2703
    Commented Nov 29, 2023 at 13:07
  • $\begingroup$ 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. $\endgroup$ Commented Nov 29, 2023 at 15:19

1 Answer 1

2
$\begingroup$

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

$\endgroup$
7
  • $\begingroup$ 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. $\endgroup$
    – je2703
    Commented Nov 29, 2023 at 14:23
  • $\begingroup$ 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. $\endgroup$
    – mmikkelsen
    Commented Nov 29, 2023 at 15:04
  • $\begingroup$ 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. $\endgroup$ Commented Nov 29, 2023 at 15:30
  • $\begingroup$ @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. $\endgroup$
    – je2703
    Commented Nov 29, 2023 at 17:04
  • 2
    $\begingroup$ 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. $\endgroup$ Commented Nov 29, 2023 at 18:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.