我正在尝试求解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 可以立即做到这一点。
怎么了?
在 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
但我们可以看到它是在舍入误差的顺序上。
您可以使用prec
参数获得更多精度数字来nsolve
:
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
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]])