Search code examples
pythonsympynonlinear-optimization

Solving system of nonlinear (and nonpolynomial) equations


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.


Solution

  • 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}