Sympy doesn't find a solution for a system of four equations when I already know the solution exists

Question

I already know that the solution to a system of four equations are these:

• 
• 
• 
• .

The SymPy code I'm using:

z1, z2, w1, w2 = symbols( 'z1, z2, w1, w2', real=True )
ex1 = 2 * (z1*w1 + z2*w2)
ex2 = w1**2 + w2**2 + 4*z1*z2*w1*w2
ex3 = 2*w1*w2 * (z1*w2 + z2*w1)
ex4 = w1**2 * w2**2
eq1 = Eq( ex1, 10 )
eq2 = Eq( ex2, 45 )
eq3 = Eq( ex3, 105 )
eq4 = Eq( ex4, 105 )
solve( [ eq1, eq2, eq3, eq4 ], [ z1, z2, w1, w2 ] )
[]                      # Note that the result is empty.

I can demonstrate the solution I have, though:

solution = { z1:0.620702965049498, z2:0.957974461859108, w1:3.38936579272158, w2:3.02326493881663 }
for i in [ ex1, ex2, ex3, ex4 ]: i.subs( solution )
9.99999999999999        # Very close to 10 in eq1
44.9999999999999        # Very close to 45 in eq2
105.000000000000        # Matches 105 in eq3
105.000000000000        # Matches 105 in eq4

I'm wondering what I may have done wrong in setting this up for SymPy. Or if there's a known solver problem? (My SymPy version is '1.10.1'.)

Note: I'm working on the expression of a homogeneous solution for a 4th order Bessel equation as used in electronics filter design, where I've factored it into two 2nd order expressions.

Answer 1

Since you are looking for numerical solutions, you probably don't need SymPy.

Here I'm going to use scipy's fsolve to solve your system of non-linear equations:

from scipy.optimize import fsolve
# convert your equations to numerical functions
eqns = [lambdify([z1, z2, w1, w2], eq.rewrite(Add)) for eq in [eq1, eq2, eq3, eq4]]

def func(p):
    # p contains the current iteration parameters: z1, z2, w1, w2
    return [eq(*p) for eq in eqns]

sol = fsolve(func, (1, 1, 1, 1))
print(sol)
# [0.62070297 0.95797446 3.38936579 3.02326494]

Note that you can also use Sympy's nsolve to solve a system of equations for numerical values, but you need to provide some good initial guess, which for your case might not be trivial:

nsolve([eq1, eq2, eq3, eq4], [z1, z2, w1, w2], [0.5, 1, 3, 3])