python - SymPy 无法求解三角方程组

Question

我试图让 SymPy 求解方程组，但它给我一个错误提示:

NotImplementedError: could not solve 3*sin(3*t0/2)*tan(t0) + 2*cos(3*t0/2) - 4

还有其他方法可以求解方程组吗:

sin(x)+(y-x)cos(x)           = 0

-1.5(y-x)sin(1.5x)+cos(1.5x) = 2

我用过:

from sympy import *
solve([sin(x)+(y-x)cos(x), -1.5(y-x)sin(1.5x)+cos(1.5x)-2], x, y)

Answer 1

SymPy 可以用这个等式做得更好，但最终它等同于某个 10 次多项式，其根只能抽象地表示。我将描述可以采取的步骤并展示 SymPy 可以走多远。这是一个应该更加自动化的半手动解决方案过程。

首先，不要在方程式中输入 1.5 或其他 float 。相反，引入一个系数 a = Rational(3, 2) 并使用它:

eq = [sin(x) + (y-x)*cos(x), -a*(y-x)*sin(a*x) + cos(a*x) - 2]

可以使用第一个方程消除变量 y:y=x-tan(x)，这对我们来说很容易看出，但 SymPy 有时会错过机会。让我们帮助它:

eq1 = eq[1].subs(y, x-tan(x))   #   3*sin(3*x/2)*tan(x)/2 + cos(3*x/2) - 2

照原样，solve 和 solveset(另一种 SymPy 求解器)由于不同参数的三角函数的这种混合而放弃方程。我们中的一些人记得在学生时代，三角函数可以表示为半参数正切的有理函数，所以让我们这样做:根据 tan 重写方程。

eq2 = eq1.rewrite(tan)   #   (-tan(3*x/4)**2 + 1)/(tan(3*x/4)**2 + 1) - 2 + 3*tan(3*x/4)*tan(x)/(tan(3*x/4)**2 + 1)

如前所述，这使论点减半了。在三角函数中有像 x/4 这样的分数是不好的。引入一个新符号 var('u')，并使 u = x/4:

eq3 = eq2.subs(x, 4*u)   #   (-tan(3*u)**2 + 1)/(tan(3*u)**2 + 1) - 2 + 3*tan(3*u)*tan(4*u)/(tan(3*u)**2 + 1)

现在我们可以使用 expand_trig 根据 tan(u) 展开所有这些切线。等式变长了:

eq4 = expand_trig(eq3)  #  (1 - (-tan(u)**3 + 3*tan(u))**2/(-3*tan(u)**2 + 1)**2)/(1 + (-tan(u)**3 + 3*tan(u))**2/(-3*tan(u)**2 + 1)**2) - 2 + 3*(-4*tan(u)**3 + 4*tan(u))*(-tan(u)**3 + 3*tan(u))/((1 + (-tan(u)**3 + 3*tan(u))**2/(-3*tan(u)**2 + 1)**2)*(-3*tan(u)**2 + 1)*(tan(u)**4 - 6*tan(u)**2 + 1))

但它也更简单，因为 tan(u) 可以被视为另一个未知数，比如 v。

eq5 = eq4.subs(tan(u), v)  #  (1 - (-v**3 + 3*v)**2/(-3*v**2 + 1)**2)/(1 + (-v**3 + 3*v)**2/(-3*v**2 + 1)**2) - 2 + 3*(-4*v**3 + 4*v)*(-v**3 + 3*v)/((1 + (-v**3 + 3*v)**2/(-3*v**2 + 1)**2)*(-3*v**2 + 1)*(v**4 - 6*v**2 + 1))

太好了，现在我们有了一个有理函数。它可以用 solveset(eq5, x) 处理。默认情况下，solveset 给出了所有复杂的解决方案，我们只需要其中的实根，所以让我们将域指定为实数:

vsol = list(solveset(eq5, v, domain=S.Reals))

这些没有代数公式，所以它们的记录有些抽象，但这些是我们可以使用的实际数字:

[CRootOf(3*v**10 + 9*v**8 - 78*v**6 + 22*v**4 - 21*v**2 + 1, 0),
 CRootOf(3*v**10 + 9*v**8 - 78*v**6 + 22*v**4 - 21*v**2 + 1, 1),
 CRootOf(3*v**10 + 9*v**8 - 78*v**6 + 22*v**4 - 21*v**2 + 1, 2),
 CRootOf(3*v**10 + 9*v**8 - 78*v**6 + 22*v**4 - 21*v**2 + 1, 3)]

例如，我们现在可以回到 x 和 y，并评估解决方案:

xsol = [4*atan(v) for v in vsol] 
ysol = [x - tan(x) for x in xsol]
numsol = [(N(x), N(y)) for x, y in zip(xsol, ysol)]

数值是

[(-4.35962510714700, -1.64344290066272),
 (-0.877886785847899, 0.326585146723377),
 (0.877886785847899, -0.326585146723377),
 (4.35962510714700, 1.64344290066272)]

当然还有无限多因为切线是周期性的。最后，让我们检查一下这些是否有效:

residuals = [[e.subs({x: xv, y: yv}) for e in eq] for xv, yv in numsol]

这些是一串 1e-15 或更小阶数，所以是的，方程式在机器精度范围内成立。

与我们从 SciPy 或其他数值求解器获得的纯数值解不同，这些可以在不重复过程的情况下以任何精度进行评估。比如第一个x-解的50位:

xsol[0].evalf(50)  #   -4.3596251071470021258397061103704574594477338857831