I'm trying to port a utility from sympy to sagemath.
Here is my sympy code:
from sympy import symbols, Eq, solve, Matrix, pprint
x, y, z = symbols('x y z')
def sympy_solve(augmented_matrix):
# Get the number of variables
num_variables = augmented_matrix.shape[1] - 1
# Generate symbols for variables
variables = symbols('x:' + str(num_variables))
# Extract coefficients and constants from the augmented matrix
coefficients = augmented_matrix[:, :-1]
constants = augmented_matrix[:, -1]
# Create equations from the coefficients and constants
equations = []
for i in range(len(constants)):
equation = Eq(sum(coefficients[i, j] * variables[j] for j in range(num_variables)), constants[i])
equations.append(equation)
# Solve the equations
solution = solve(equations, variables, dict=True)
return solution
When I run it, I get the result I expect:
A = Matrix([
[1,1,1,-1,1],
[0,1,-1,1,-1],
[3,0,6,-6,6],
[0,-1,1,-1,1]
])
print(sympy_solve(A))
# [{
# x0: -2*x2 + 2*x3 + 2,
# x1: x2 - x3 - 1
# }]
My current solution in sagemath looks like this:
def sagemath_solve(augmented_matrix):
A = augmented_matrix[:, :-1]
Y = augmented_matrix[:, -1]
m, n = A.dimensions()
p, q = Y.dimensions()
if m!=p:
raise RuntimeError("The matrices have different numbers of rows")
X = vector([var("x_{}".format(i)) for i in [1..n]])
sols = []
for j in range(q):
system = [A[i]*X==Y[i,j] for i in range(m)]
sols += solve(system, *X)
return sols
However, it is returning rationals:
A = matrix([
[1,1,1,-1,1],
[0,1,-1,1,-1],
[3,0,6,-6,6],
[0,-1,1,-1,1]
])
print(sagemath_solve(A))
# [[
# x_1 == 2*r1 - 2*r2 + 2,
# x_2 == -r1 + r2 - 1,
# x_3 == r2,
# x_4 == r1
# ]]
Any thoughts on how to fix this?