2

I am trying to solve a system of multivariate equations, which are the result of some Java code. Neither the form, nor the number of variables is known before runtime. An example would be

(I) (e-a*d*e-b*d*e+2*b*d*f+2*b*d*e*g)/(-1+a*d+b*d)+f == 0

(II) e*g+((f+e*g)*a*d)/(-1+a*d+b*d)==0

(III) -e*h+((-f-e*g)*d)/(-1+a*d+b*d)==0

(IV) -e*j+((-f-e*g)*c)/(-1+a*d+b*d)==0

I tried using Symja, which simply returns the input, and SymPy, which throws an Error

ZeroDivisionError: polynomial division

The variables are all from the interval [0,1], and I need all solutions. Mathematica is able to solve this, but as it is commercial software I unfortunately cannot use it in this project.

I would be grateful for any recommendations on which software to use. I really would have liked SymPy to work, and I don't understand why it throws this error, ideas are appreciated. Below a MWE to the SymPy Error:

from sympy.solvers import solve
from sympy.abc import a,b,c,d,e,f,g,h,j

lst = a,b,c,d,e,f,g,h,j
sys = [(e-a*d*e-b*d*e+2*b*d*f+2*b*d*e*g)/(-1+a*d+b*d)+f,e*g+((f+e*g)*a*d)/(-1+a*d+b*d),-e*h+((-f-e*g)*d)/(-1+a*d+b*d),-e*j+((-f-e*g)*c)/(-1+a*d+b*d)]

solution = solve(sys, lst)
print solution

The Mathematica version is:

eqn = {(e - a*d*e - b*d*e + 2*b*d*f + 2*b*d*e*g)/(-1 + a*d + b*d) + f == 0, e*g + ((f + e*g)*a*d)/(-1 + a*d + b*d) == 0, -e*h + ((-f - e*g)*d)/(-1 + a*d + b*d) == 0, -e*j + ((-f - e*g)*c)/(-1 + a*d + b*d) == 0};
Simplify[Solve[eqn, {a, b, c, d, e, f, g, h, j}]]

Output:

{{e -> 0, f -> 0},
{c -> (1 - 2 a d - 3 b d) j, f -> ((-1 + 2 a d + b d) e)/(-1 + 2 a d + 3 b d), g -> (a d)/(1 - 2 a d - 3 b d), h -> d/(1 - 2 a d - 3 b d)},
{a -> 0, c -> j - 3 b d j, f -> ((-1 + b d) e)/(-1 + 3 b d), g -> 0, h -> d/(1 - 3 b d)},
{a -> (1 - b d)/(2 d), c -> -2 b d j, f -> 0, g -> 1/4 - 1/(4 b d), h -> -(1/(2 b))}}
berndibus
  • 95
  • 10

3 Answers3

1

Note that you have 9 variables and 4 equations. Hence, you can eliminate four variables - you need to tell Sympy which.

The following code first multiplies the denominator out and the eliminates a,b,c,d:

import sympy as sy
from IPython.display import display  # for pretty printing

# sy.init_printing()  # LaTeX-like pretty printing for IPython

a, b, c, d, e, f, g, h, j = sy.symbols("a, b, c, d, e, f, g, h, j", real=True)
lst = a, b, c, d, e, f, g, h, j
sys0 = sy.Matrix([(e-a*d*e-b*d*e+2*b*d*f+2*b*d*e*g)/(-1+a*d+b*d)+f,
                  e*g+((f+e*g)*a*d)/(-1+a*d+b*d),
                  -e*h+((-f-e*g)*d)/(-1+a*d+b*d),
                  -e*j+((-f-e*g)*c)/(-1+a*d+b*d)])

# Denominator can be factored out:
den = a*d + b*d - 1
sys1 = sy.simplify(sys0*den).expand()
print("Factored out denominator:")
display(sys1)

# Elimnate four variables:
sol1 = sy.solve(sys1, a, b, c, d, dict=True)
print("Solutions:")
display(sol1)
print("Substituting back into the equation gives obviously 0, i.e.:")
display(sy.simplify(sys1.subs(sol1[0])).T)

print("The denominator != 0 results in:")
den1 = sy.simplify(den.subs(sol1[0]))
display(sy.solve(den1))

which generates the following output:

Factored out denominator:
Matrix([
[-a*d*e + a*d*f + 2*b*d*e*g - b*d*e + 3*b*d*f + e - f],
[                   2*a*d*e*g + a*d*f + b*d*e*g - e*g],
[              -a*d*e*h - b*d*e*h - d*e*g - d*f + e*h],
[              -a*d*e*j - b*d*e*j - c*e*g - c*f + e*j]])
Solutions:
[{d: 2*e*h/(4*e*g - e + 3*f),
  c: 2*e*j/(4*e*g - e + 3*f),
  a: g/h,
  b: (-e + f)/(2*e*h)}]
Substituting back into the equation gives obviously 0, i.e.:
Matrix([[0, 0, 0, 0]])
The denominator != 0 results in:
[{e: -f/g}]

So the resulting identities are:

-f/g != e  # denominator - only false, if f,e=0 since f,g,e>=0
a = g/h
b = (-e + f)/(2*e*h)
c = 2*e*j/(4*e*g - e + 3*f)
d = 2*e*h/(4*e*g - e + 3*f)

As far as I'm aware, Sympy cannot handle multivariate inequalities yet (though I would like to be proven wrong). But the results are simple enough to do by hand:

 0 <= g <= h
 0 <= f-e <= 2*e*h
 0 <= 2*e*j <= 4*e*g - e + 3*f
 0 <= 2*e*h <= 4*e*g - e + 3*f
 0 <= e,f,g,h,j <= 1

The case f,e=0 is valid as well. It can be verified by checking that sys0.subs(e,0).diff(f) does not depend on f.

Dietrich
  • 5,241
  • 3
  • 24
  • 36
  • 1
    very nice approach. – Stelios Jan 05 '17 at 17:13
  • That is indeed a nice approach for this particular system. Unfortunately, the structure of the equations are not known before runtime. That means that this case of all equations having a common denominator is an artifact, and may be different in other cases. Thus I do not see how this solution can be extended to work programmatically. – berndibus Jan 06 '17 at 00:39
  • @berndibus If you know nothing or little about the structure of your equations, in my experience such an approach is problematic (especially for the corner cases). In Mathematica (and also Sympy), I routinely wait half an hour and more to see if `Solve[]` comes up with a valid solution or not . Note that Mathematica does not give *all* solutions in your example - `Reduce[]` would be more appropriate. If you're sure, you're equations are well behaved, your 80% solution might be to eliminate the first four variables as I done in my solution. – Dietrich Jan 06 '17 at 09:03
0

In case this question gets closed, here is an attempt at the Mathematica version. However it currently lacks conditions restricting the variables to the interval [0,1].

Solve[{
  (e - a d e - b d e + 2 b d f + 2 b d e g)/(-1 + a d + b d) + f == 0,
   e g + ((f + e g) a d)/(-1 + a d + b d) == 0,
  -e h + ((-f - e g) d)/(-1 + a d + b d) == 0,
  -e j + ((-f - e g) c)/(-1 + a d + b d) == 0
  },
 {a, b, c, d, e, f, g, h, j}]

enter image description here

Chris Degnen
  • 8,443
  • 2
  • 23
  • 40
  • I have added the Mathematica version I have to my question, but I cannot use Mathematica for my problem as it is commercial software. – berndibus Jan 04 '17 at 14:57
  • You have not added any conditions that would limit the variables to the interval [0,1] so guess that is not primarily relevant to what you require from SymPy (or other software). – Chris Degnen Jan 04 '17 at 15:20
  • It would be nice, but I thought it would be possible to restrict that in a second step, the main task is the solving. If it is possible to do that directly, it would be much better. – berndibus Jan 04 '17 at 17:12
0

Ok, I was able to solve the problem myself using SAGE, with gives the same output as Mathematica.

sage: a,b,c,d,e,f,g,h,j = var('a,b,c,d,e,f,g,h,j')
sage: qe = [(e-a*d*e-b*d*e+2*b*d*f+2*b*d*e*g)+f*(-1+a*d+b*d),e*g*(-1+a*d+b*d)+((f+e*g)*a*d),-e*h*(-1+a*d+b*d)+((-f-e*g)*d),-e*j*(-1+a*d+b*d)+((-f-e*g)*c)]
sage: print(solve(qe,a,b,c,d,e,f,g,h,j, solution_dict=True))

Gives the output

[{g: r5, j: r7, b: r2, d: r4, e: 0, h: r6, c: r3, f: 0, a: r1},
{g: r12, j: r14, b: r8, d: r10, e: r11, h: r13, c: r9, f: -r11*r12, a: -(r10*r8 - 1)/r10},
{g: r15*r18, j: r19, b: -1/6*(4*r15*r16*r18 - r16 + 3*r17)*(4*r15*r16*r18^2/(4*r15*r16*r18 - r16 + 3*r17) + 2*r16*r18/(4*r15*r16*r18 - r16 + 3*r17) - r18)/(r16*r18^2), d: 2*r16*r18/(4*r15*r16*r18 - r16 + 3*r17), e: r16, h: r18, c: 2*r16*r19/(4*r15*r16*r18 - r16 + 3*r17), f: r17, a: r15},
{g: r20*r23, j: 0, b: -1/6*(4*r20*r21*r23 - r21 + 3*r22)*(4*r20*r21*r23^2/(4*r20*r21*r23 - r21 + 3*r22) + 2*r21*r23/(4*r20*r21*r23 - r21 + 3*r22) - r23)/(r21*r23^2), d: 2*r21*r23/(4*r20*r21*r23 - r21 + 3*r22), e: r21, h: r23, c: 0, f: r22, a: r20},
{g: -2*r24*r27 - 1, j: 0, b: r24, d: r25, e: r26, h: r27, c: 0, f: 2*r24*r26*r27 + r26, a: -(r24*r25 - 1)/r25},
{g: 0, j: r30, b: r29, d: 0, e: r31, h: 0, c: r30, f: r31, a: r28},
{g: 0, j: 0, b: r33, d: 0, e: r34, h: 0, c: 0, f: r34, a: r32},
{g: -1, j: 0, b: 0, d: r35, e: r36, h: r37, c: 0, f: r36, a: 1/r35},
{g: r40, j: 0, b: 0, d: r38, e: r39, h: 2*r38*r40 + r38, c: 0, f: r39, a: r40/(2*r38*r40 + r38)},
{g: -1/2, j: 0, b: -1/2*(r41 - r42)/(r41*r43), d: -2/3*r41*r43/(r41 - r42), e: r41, h: r43, c: 0, f: r42, a: -1/2/r43}]
berndibus
  • 95
  • 10