-1

I'm trying to solve a system of equations using sympy.

from sympy import *
def get_angles(a, b, c, d):
    theta, phi, lamb = symbols('\\theta \\phi \\lambda', real=True)
    a_eq = Eq(cos(theta / 2), a)
    b_eq = Eq(exp(I * phi) * sin(theta / 2), b)
    c_eq = Eq(-exp(I * lamb) * sin(theta / 2), c)
    d_eq = Eq(exp(I * (phi + lamb)) * cos(theta / 2), d)
    # theta_constr1 = Eq(theta >= 0)
    # theta_constr2 = Eq(theta <= pi)
    # phi_constr1 = Eq(phi >= 0)
    # phi_constr2 = Eq(phi < 2 * pi)
    res = solve([
        a_eq, b_eq, c_eq, d_eq,
        #theta_constr1, theta_constr2, phi_constr1, phi_constr2,
    ],
        theta,
        phi,
        lamb,
        check=False,
        dict=True)
    return res

The function returns the right results as it is, but it doesn't work if I try to put the constraint on the angles inside the system of equations (the commented parts). Is there any way to have them?

At the moment, I'm using a simple solution to overcome this limitation: I pass the result of the previous function to the following one to filter out unwanted results

def _final_constraint(result):
    res = []
    for sol in result:
        to_add = True
        for k, v in sol.items():
            if str(k) == '\\theta' and (v < 0 or v > pi):
                to_add = False
                break
            elif str(k) == '\\phi' and (v < 0 or v >= 2 * pi):
                to_add = False
                break
        if to_add:
            res.append(simplify(sol))
    return res
tigerjack
  • 1,158
  • 3
  • 21
  • 39

1 Answers1

1

Eq stands for Equality and is not an equation (though there is discussion of adding such an object to SymPy). So your uncommented Eq get interpreted as you intended, but the commented ones don't. You could try replacing theta_constr1 = Eq(theta >= 0) with theta_constr1 = theta >= 0 but then you will run into problems with the the inequality solver -- it complains about there being more than one symbol of interest amongst the inequalities. So what about re-writing the inequalities like x >= 0 as Eq(x + eps, 0) where eps is a nonnegative Symbol:

def get_angles(a, b, c, d):
    theta, phi, lamb = symbols('\\theta \\phi \\lambda', real=True)
    eps = Symbol('eps', nonnegative=True)
    a_eq = Eq(cos(theta / 2), a)
    b_eq = Eq(exp(I * phi) * sin(theta / 2), b)
    c_eq = Eq(-exp(I * lamb) * sin(theta / 2), c)
    d_eq = Eq(exp(I * (phi + lamb)) * cos(theta / 2), d)
    theta_constr1 = theta + eps
    theta_constr2 = pi - theta + eps
    phi_constr1 = phi  + eps
    phi_constr2 = 2 * pi - phi + eps
    res = solve([
        a_eq, b_eq, c_eq, d_eq,
        theta_constr1, theta_constr2, phi_constr1, phi_constr2,
        ],
        theta,
        phi,
        lamb,
        check=False,
        set=True)
    return res

>>> print(filldedent(get_angles(11,2,3,4)))

([\lambda, \phi, \theta], {(pi, 0, 4*pi - 2*acos(11)), (0, pi,
2*acos(11)), (0, pi, 4*pi - 2*acos(11)), (0, 0, 4*pi - 2*acos(11)),
(0, 0, 2*acos(11)), (pi, pi, 4*pi - 2*acos(11)), (pi, pi, 2*acos(11)),
(pi, 0, 2*acos(11))})

You would have to work out which side of theta and phi -- if any -- satisfied your equations.

smichr
  • 16,948
  • 2
  • 27
  • 34
  • However, your solution doesn't seem to work for, say, `get_angles(.5*sqrt(2), -0.5*sqrt(2)*I, -0.5*sqrt(2)*I, 0.5*sqrt(2)`. One of the result is `{\theta: 10.9955742875643, \lambda: 1.57079632679490, \phi: 1.57079632679490}`. – tigerjack May 28 '19 at 19:44
  • If you take out `check=False` the solutions will be checked and you will see that in this case there is no solution (if I am reading the results correctly). Either an Eq is false or eps would have to be negative when we required that it not be so. – smichr May 28 '19 at 21:25
  • Well, there is the solution `{\theta: pi/2, \lambda: pi/2, \phi: 3*pi/2}`, which is never returned if I'm using those constraints. – tigerjack May 29 '19 at 07:26
  • Hmm. Looks like you should just leave the auxiliary conditions off and filter your solutions to the other equations with those conditions. Answer has been updated to show that. – smichr May 29 '19 at 17:10