Solving system of nonlinear (and nonpolynomial) equations

Question

I would like to use python/sympy to solve simple systems of equations coming from impedance calculations in electronics.

In such calculations, due to the "parallel impedance" formula, one often has to deal with expressions of the form:

par(x,y) := (x*y)/(x+y)

Now I have tried with the following code:

from sympy import *

def par(var1,var2):
    return (var1 * var2)/(var1+var2)

A = Symbol('A')
B = Symbol('B')
C = Symbol('C')
D = Symbol('D')
E = Symbol('E')

eq1 = A + par(B+50 , (C+  par(D,  (E+50))  ))  - 50 
eq2 = B + par(A+50 , (C+  par(D ,  (E+50)) ))  - 50
eq3 = E + par(D+50,  (C+ par(A+50, B+50)) ) - 50

thus defining a system of three equations in five variables {A,B,C,D,E}, but then running

solve([eq1,eq2,eq3], A, B,C,D,E)

the computations just does not terminate.

Do you have any suggestions on how I could approach these type of equations? Basically polynomials with division by polynomials, with solutions in the complex numbers.

Is this an Underdetermined system? It may be there are infinite solutions. — Perplexabot, Dec 16 '22 at 21:14
yes it is. Is there a way to ask the procedure to return with just the first solution found? — makkiato, Dec 16 '22 at 21:22
If you don't care about representing the infinite solution set then choose arbitrary values for D and E and then solve for A, B, and C. — Oscar Benjamin, Dec 16 '22 at 21:44

JohanC · Answer 1 · 2022-12-17T00:05:07.150

You could also try the Z3 library:

from z3 import Reals, Solver, sat, set_option

def par(var1, var2):
    return (var1 * var2) / (var1 + var2)

vars = Reals('A B C D E')
A, B, C, D, E = vars

set_option(rational_to_decimal=True, precision=30)
s = Solver()
s.add(A + par(B + 50, (C + par(D, (E + 50)))) == 50)
s.add(B + par(A + 50, (C + par(D, (E + 50)))) == 50)
s.add(E + par(D + 50, (C + par(A + 50, B + 50))) == 50)
if s.check() == sat:
    m = s.model()
    print(m)

I get following output:

[D = 0.502507819412956050111927878839?,
 C = 0.125,
 E = 50.188198722936229689352204927638?,
 B = -50.625,
 A = -50.625]

The question mark at the end of D's and E's values means they have been approximated.

If you then try the values for A, B and C into the original code, sympy gives two exact expressions:

[{A: -405/8,
  B: -405/8,
  C: 1/8,
  D: -31857/1292 + 5*sqrt(676050154)/5168,
  E: 443/1304 + 5*sqrt(676050154)/2608},
 {A: -405/8,
  B: -405/8,
  C: 1/8,
  D: -5*sqrt(676050154)/5168 - 31857/1292,
  E: 443/1304 - 5*sqrt(676050154)/2608}]

smichr · Accepted Answer · 2022-12-17T03:33:38.410

Taking the suggestion of Oscar and focussing on A, B and C you can get a solution for them in terms of D and E:

>>> solve((eq1, eq2, eq3), A, B, C)[0] # one solution
(
50*(D + E)*(D + E + 50)/(3*D**2 - 2*D*E + 150*D - 3*E**2 - 50*E + 5000), 
50*(D + E)*(D + E + 50)/(3*D**2 - 2*D*E + 150*D - 3*E**2 - 50*E + 5000), 
(-3*D**3*E + 50*D**3 + 2*D**2*E**2 - 500*D**2*E + 10000*D**2 + 3*D*E**3 + 50*D*E**2 - 25000*D*E + 500000*D + 200*E**3 - 5000*E**2 - 500000*E + 12500000)/(3*D**3 + D**2*E + 150*D**2 - 5*D*E**2 + 100*D*E + 5000*D - 3*E**3 - 50*E**2 + 5000*E))

Notice that the solution for A and B is the same (consistent with the symmetry of the first two equations wrt A and B).

>>> sol = Dict(*zip((A,B,C),_))
>>> sol[A] = sol[B]
True

In this form, you can directly substitute values for D and E:

>>> sol.subs({D:1, E:S.Half})
{A: 1030/1367, B: 1030/1367, C: 2265971/1367}

You can also see what relationships between D and E are forbidden by solving for when any of the denominators are 0:

>>> from sympy.solvers.solvers import denoms
>>> set([j.simplify() for i in denoms(sol) for j in solve(i,D)])
{-E, E/3 - sqrt(10*E**2 - 9375)/3 - 25, E/3 + sqrt(10*E**2 - 9375)/3 - 25}

Solving system of nonlinear (and nonpolynomial) equations

2 Answers2