How to solve non-linear system of trigonometric equations in Python (which MATLAB can solve easily)

Question

I am trying to solve a system of non-linear trigonometric equations in Python. I tried the following:

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 takes forever to attempt to solve the equations. However, MATLAB does it in an instant.

What's the problem?

Then Python should also give a non-analytical solution. But it doesn't give anything. — Swapnil Saha, Dec 13 '19 at 19:32
In MATLAB, you can combine symbolic mathematics with numerical solvers. Python doesn't seem to support that — Swapnil Saha, Dec 13 '19 at 19:39
Could you provide a function that will provide a numerical solution? — Swapnil Saha, Dec 13 '19 at 19:40
Python is a general purpose programming language. What you do with it depends on the imported modules. — , Dec 13 '19 at 19:40

score 6 · Accepted Answer · answered Dec 13 '19 at 20:37

In SymPy if you want numerical solutions you should use nsolve:

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

Here I've used an initial guess of [1, 1, 1]. I'm sure you can find more solutions if you use other initial guesses (the system has an infinite number of solutions).

Note that if you substitute these approximate solutions into the equations you will get False. That's because the lhs and the rhs as approximate numbers are unequal:

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

Since that isn't equal to the rhs (which is zero) substituting into the Eq will give False but we can see that it is on the order of rounding error.

You can get more digits of accuracy using the prec argument to 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

score 2 · Answer 2 · answered Dec 13 '19 at 20:33

2

Sympy can also search for numerical solutions, so you can keep the format the equations are is. Note that nsolve internally uses the multiprecission library mpmath and needs a set of initial values.

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)))

Output:

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

answered Dec 13 '19 at 20:33

JohanC

71,591
8
33
66

Thanks. This works. The only problem here is the initial guess. It, especially the guess for f (the frequency of the signal). I guess i will need to find some good initial guesses for it to come close to what I want. – Swapnil Saha Dec 13 '19 at 21:22
1

You can solve the first two equations for f and phi with `solve([eq1, eq2], [f, phi], dict=True)`. That gives 4 analytic solutions for f and phi in terms of a. Then if you substitute one of those into the third equation you can solve numerically for `a` with `nsolve`. At that point there is a unique solution for `a` so an initial guess of 1 will probably always work. – Oscar Benjamin Dec 13 '19 at 22:32
@SwapnilSaha: Please change the accepted solution to OscarBenjamin's, as it is more thorough than mine. – JohanC Dec 13 '19 at 22:47

How to solve non-linear system of trigonometric equations in Python (which MATLAB can solve easily)

2 Answers2