Welcome to Software Development on Codidact!
Will you help us build our independent community of developers helping developers? We're small and trying to grow. We welcome questions about all aspects of software development, from design to code to QA and more. Got questions? Got answers? Got code you'd like someone to review? Please join us.
How to apply a unit constraint in SymPy
The solution to this is probably easy, but I haven't been able to find it. Using SymPy, I am trying to solve this equation: $$ {x_1}^2 + {x_2}^2 + {x_3}^2 + 3 = a $$ with this constraint: $$ {x_1}^2 + {x_2}^2 + {x_3}^2 = 1 $$ It is obvious that $a=4$. So, I wrote this code:
import sympy as sp
a, x1, x2, x3 = sp.symbols('a x1 x2 x3')
eq1 = x1**2 + x2**2 + x3**2 + 3 - a
eq2 = x1**2 + x2**2 + x3**2 - 1
sol = sp.solve([eq1, eq2], a)
print(sol)
Which returns:
{a: x1**2 + x2**2 + x3**2 + 3}
The unit-norm constraint on x is not applied.
How can I get SymPy to apply the unit-norm constraint, which results in $a=4$?
3 answers
From the documentation for solve:
If any equation does not depend on the symbol(s) given, it will be eliminated from the equation set and an answer may be given implicitly in terms of variables that were not of interest:
>>> solve([x - y, y - 3], x) {x: y}
That's what you're encountering here, as your eq2 doesn't depend on a.
Simplest thing is to omit the symbol altogether:
>>> sp.solve([eq1, eq2])
[{a: 4, x1: -sqrt(-x2**2 - x3**2 + 1)}, {a: 4, x1: sqrt(-x2**2 - x3**2 + 1)}]
Unfortunately, this requires a bit of extra interpretation to see that the value found for a is the same in both solutions. I think that's just something to live with.
0 comment threads
To expand on r's answer, the "bit of extra interpretation" might look like this:
r = sp.solve([eq1, eq2])
# Solution set
a_solutions = {d[a] for d in r}
# Assert one solution
assert len(a_solutions) == 1
# Extract it
a_solved = next(iter(a_solutions))
0 comment threads
Are the coefficients of eq1 always 1? If so, you can simply subtract one equation from the other.
>>> sp.solve(eq1-eq2, a)
[4]
P.S. I don't have much experience with SymPy, so LMK if I missed some subtlety.
1 comment thread
0 comment threads