如何在 Python 中求解非线性三角方程组（MATLAB 可以轻松求解）答案

【问题标题】：How to solve non-linear system of trigonometric equations in Python (which MATLAB can solve easily)如何在 Python 中求解非线性三角方程组（MATLAB 可以轻松求解）
【发布时间】：2019-12-13 19:14:06
【问题描述】：

我正在尝试在 Python 中求解非线性三角方程组。我尝试了以下方法：

from sympy import symbols,solve,sin,cos,pi, Eq

measurements = [(5.71403,0.347064), (4.28889, -0.396854), (5.78091, -7.29133e-05), 
(2.06098, 0.380579), (8.13321, 0.272391), (8.23589, -0.304111), (6.53473, 0.265354), (1.6023, 
0.131908)]

f, a, phi = symbols('f a phi')
eq1 = Eq(a*sin((2.0*pi*f*measurements[0][0])+phi) - measurements[0][1])
eq2 = Eq(a*sin((2.0*pi*f*measurements[4][0])+phi) - measurements[4][1])
eq3 = Eq(a*sin((2.0*pi*f*measurements[6][0])+phi) - measurements[6][1])
solve((eq1,eq2,eq3), (a, f, phi))

Python 需要很长时间才能尝试求解方程。然而，MATLAB 瞬间完成。

有什么问题？

【问题讨论】：

我怀疑是否有解析解。
当时MATLAB是怎么解决的？
显然，不是分析性的。
那么Python也应该给出一个非解析解。但它什么也没给。
Python 是一种通用编程语言。你用它做什么取决于导入的模块。

标签： python numpy scipy sympy

【解决方案1】：

在 SymPy 中，如果你想要数值解，你应该使用 nsolve:

In [97]: nsolve((eq1,eq2,eq3), (a, f, phi), [1, 1, 1])                                                                            
Out[97]: 
⎡-0.5538674055548 ⎤
⎢                 ⎥
⎢0.837453526933376⎥
⎢                 ⎥
⎣6.95538865037068 ⎦

这里我使用了[1, 1, 1] 的初始猜测。如果您使用其他初始猜测，我相信您可以找到更多解决方案（系统有无限数量的解决方案）。

请注意，如果您将这些近似解代入方程，您将得到 False。那是因为作为近似数的 lhs 和 rhs 是不相等的：

In [101]: eq1                                                                                                                     
Out[101]: a⋅sin(11.42806⋅π⋅f + φ) - 0.347064 = 0

In [102]: (sol,) = nsolve((eq1,eq2,eq3), (a, f, phi), [1, 1, 1], dict=True)                                                       

In [103]: sol                                                                                                                     
Out[103]: {a: -0.5538674055548, f: 0.837453526933376, φ: 6.95538865037068}

In [104]: eq1.subs(sol)                                                                                                           
Out[104]: False

In [105]: eq1.lhs.subs(sol)                                                                                                       
Out[105]: -0.347064 - 0.5538674055548⋅sin(6.95538865037068 + 9.57046915300624⋅π)

In [106]: eq1.lhs.subs(sol).evalf()                                                                                               
Out[106]: -1.29025679909939e-15

由于不等于 rhs（为零）代入 Eq 将得到 False，但我们可以看到它是舍入误差的顺序。

您可以使用nsolve 的prec 参数获得更多位数的准确度：

In [107]: (sol,) = nsolve((eq1,eq2,eq3), (a, f, phi), [1, 1, 1], dict=True, prec=50)                                              

In [108]: sol                                                                                                                     
Out[108]: 
{a: -0.55386740555480009188439615822304411607289430639164, f: 0.83745352693337644862065403386504543698722276260565, φ: 6.9553886
503706758809942541544797040214354242211993}

In [109]: eq1.lhs.subs(sol).evalf()                                                                                               
Out[109]: -3.27785083138700e-51

【讨论】：

【解决方案2】：

Sympy 还可以搜索数值解，因此您可以保持方程的格式。请注意，nsolve 内部使用了多精度库 mpmath，并且需要一组初始值。

from sympy import symbols, sin, pi, Eq, nsolve

measurements = [(5.71403,0.347064), (4.28889, -0.396854), (5.78091, -7.29133e-05),
(2.06098, 0.380579), (8.13321, 0.272391), (8.23589, -0.304111), (6.53473, 0.265354), (1.6023,
0.131908)]

f, a, phi = symbols('f a phi')
eq1 = Eq(a*sin((2.0*pi*f*measurements[0][0])+phi) - measurements[0][1])
eq2 = Eq(a*sin((2.0*pi*f*measurements[4][0])+phi) - measurements[4][1])
eq3 = Eq(a*sin((2.0*pi*f*measurements[6][0])+phi) - measurements[6][1])
print(nsolve((eq1,eq2,eq3), (a, f, phi), (1, 1, 0)))

输出：

Matrix([[-0.677229584607299], [1.64528629772987], [-23.9739925277907]])

【讨论】：

谢谢。这行得通。这里唯一的问题是最初的猜测。它，尤其是对 f（信号的频率）的猜测。我想我需要找到一些好的初步猜测才能接近我想要的。
您可以使用solve([eq1, eq2], [f, phi], dict=True) 求解 f 和 phi 的前两个方程。就 a 而言，这给出了 f 和 phi 的 4 个解析解。然后，如果您将其中一个代入第三个等式，您可以用nsolve 对a 进行数值求解。那时，a 有一个独特的解决方案，因此初始猜测 1 可能总是有效。
@SwapnilSaha：请将接受的解决方案更改为 OscarBenjamin 的，因为它比我的更彻底。