Sympy是一个完全基于Python的CAS，所以我不相信像求解器这样的复杂计算可以轻松地移植到GPU上。
然而，通过查看您的示例，似乎您正在尝试获得数值解。因此，solve不是正确的工具。nsolve（数值解）可能更适合这项工作。通常，您需要提供适当的初始猜测。例如：

l_bound = 4
u_bound = 13
expr1 = y/(x-1) - 3 * y/sqrt(((x-1)**2)*(x-2))-l_bound
expr2 = y/(x-1) + 3 * y/sqrt(((x-1)**2) * (x-2))-u_bound
nsolve([expr1, expr2], [x, y], [10, 10])
# out: (34.1111111111111, 281.444444444444)

你会注意到，对于这个特殊的例子，nsolve比solve快得多。然而，对于你的循环的每次迭代，nsolve将使用lambdify将符号方程组转换为数值函数组：这个操作是昂贵的。如果你的方程的形式没有随着迭代而改变，那么我们可以在开始时将符号方程系统转换为数值函数系统一次，然后在循环中使用更新的值重新使用它。例如：

from mpmath import findroot

l_bound, u_bound = symbols("l u")
expr1 = y/(x-1) - 3 * y/sqrt(((x-1)**2)*(x-2))-l_bound
expr2 = y/(x-1) + 3 * y/sqrt(((x-1)**2) * (x-2))-u_bound

system = Matrix([expr1, expr2]).T
J = system.jacobian([x, y])
# convert symbolic objects to numerical functions
system = lambdify([x, y, l_bound, u_bound], system.T, modules="mpmath")
J = lambdify([x, y, l_bound, u_bound], J, modules="mpmath")

l_bounds, u_bounds = [4], [13]
IG = [10, 10]
for l, u in zip(l_bounds, u_bounds):
    # provide new values of lower and upper bounds with each iteration
    sol = findroot(lambda x, y: system(x, y, l, u), IG, J=lambda x, y: J(x, y, l, u))
    # continue your task

备注：

1. sympy的nsolve将完全执行我们在前面代码块中所做的操作：将方程组转换为符号矩阵，计算符号雅可比矩阵，将系统和雅可比矩阵转换为数值函数，并使用mpmath的findroot来找到解决方案。
   1.如果你不使用mpmath的findroot，而是使用Scipy的一些其他的求根算法，也许可以进一步提高性能。