Trigonometric simultaneous equation

Question

I can't solve this trigonometric simultaneous equation.

(1) cos(C)=-sin(B)*sin(D)*cos(E)+sin(D)*sin(E)*cos(B)

(2) -sin(A)*sin(B)*sin(D)*sin(E)-sin(A)*sin(D)*cos(B)*cos(E)-cos(A)*cos(D)=0

I'd like to get sin(D), sin(E) only using the angle A, B, C which are constant numbers.

I tried the code below and got the results in the figure which contains cos (D) and cos(E). Variables D and E need to be eliminated. How should I do for this?

import sympy as sp
from sympy import sin, cos
sp.init_printing()
sp.var('A,B,C,D,E')
eq1=sp.Eq(cos(C),-sin(B)*sin(D)*cos(E)+sin(D)*sin(E)*cos(B))
eq2=sp.Eq(0,-sin(A)*sin(B)*sin(D)*sin(E)-sin(A)*sin(D)*cos(B)*cos(E)-cos(A)*cos(D))
sp.solve ([eq1, eq2], [sin(D), sin(E)]) 

Answer 1

I think that the result you are looking for here is more complicated than you might expect. Your equations are polynomial in sin(A), cos(A) etc. We also have polynomial relations like sin(A)**2 + cos(A)**2 - 1 = 0. Putting these together we can systematically eliminate variables using the method of resultants. Here I've replaced sin(A) and cos(A) with symbols sA and cA so that we can more clearly see that we are working with polynomials:

from sympy import *

sA, cA, sB, cB, sC, cC, sD, cD, sE, cE = symbols('sA, cA, sB, cB, sC, cC, sD, cD, sE, cE')

eq1 = Eq(cC, -sB*sD*cE + sD*sE*cB)
eq2 = Eq(0, -sA*sB*sD*sE - sA*sD*cB*cE - cA*cD)

p1 = eq1.rewrite(Add)
p2 = eq2.rewrite(Add)

# Eliminate cD
p3 = sD**2 + cD**2 - 1
p1_1 = resultant(p1, p3, cD)
p2_1 = resultant(p2, p3, cD)

# Eliminate cE
p4 = sE**2 + cE**2 - 1
p1_2 = resultant(p1_1, p4, cE)
p2_2 = resultant(p2_1, p4, cE)

# Now eliminate sD to get an univariate polynomial for sE
p_final = resultant(p1_2, p2_2, sD)

This now gives us a polynomial equation for sE with coefficients that include cA,cB,cC,sA,sB i.e. we have eliminated cos(E), sin(D) and cos(D) from the system reducing it down to a single polynomial equation for the single unknown sin(E). The roots of this polynomial will give us the values of sin(E) that would be the solutions of the system of equations that you showed (taking into account the important relationships between sin(D) and cos(D) and so on). The polynomial for sE is of degree 16 and the expressions for the coefficients are actually too complicated to be included in a SO post which is limited to 30000 characters. Instead I will just show the leading coefficient (the coefficient of sE**16):

In [47]: p = Poly(p_final, sE)

In [48]: str(p.LC())
Out[48]: 'cA**16*cB**16 + 8*cA**16*cB**14*sB**2 + 28*cA**16*cB**12*sB**4 + 56*cA**16*cB**10*sB**6 + 70*cA**16*cB**8*sB**8 + 56*cA**16*cB**6*sB**10 + 28*cA**16*cB**4*sB**12 + 8*cA**16*cB**2*sB**14 + cA**16*sB**16 + 8*cA**14*cB**16*cC**2*sA**2 + 32*cA**14*cB**14*cC**2*sA**2*sB**2 + 32*cA**14*cB**12*cC**2*sA**2*sB**4 - 32*cA**14*cB**10*cC**2*sA**2*sB**6 - 80*cA**14*cB**8*cC**2*sA**2*sB**8 - 32*cA**14*cB**6*cC**2*sA**2*sB**10 + 32*cA**14*cB**4*cC**2*sA**2*sB**12 + 32*cA**14*cB**2*cC**2*sA**2*sB**14 + 8*cA**14*cC**2*sA**2*sB**16 + 28*cA**12*cB**16*cC**4*sA**4 + 32*cA**12*cB**14*cC**4*sA**4*sB**2 - 112*cA**12*cB**12*cC**4*sA**4*sB**4 - 288*cA**12*cB**10*cC**4*sA**4*sB**6 - 344*cA**12*cB**8*cC**4*sA**4*sB**8 - 288*cA**12*cB**6*cC**4*sA**4*sB**10 - 112*cA**12*cB**4*cC**4*sA**4*sB**12 + 32*cA**12*cB**2*cC**4*sA**4*sB**14 + 28*cA**12*cC**4*sA**4*sB**16 + 56*cA**10*cB**16*cC**6*sA**6 - 32*cA**10*cB**14*cC**6*sA**6*sB**2 - 288*cA**10*cB**12*cC**6*sA**6*sB**4 + 32*cA**10*cB**10*cC**6*sA**6*sB**6 + 464*cA**10*cB**8*cC**6*sA**6*sB**8 + 32*cA**10*cB**6*cC**6*sA**6*sB**10 - 288*cA**10*cB**4*cC**6*sA**6*sB**12 - 32*cA**10*cB**2*cC**6*sA**6*sB**14 + 56*cA**10*cC**6*sA**6*sB**16 + 70*cA**8*cB**16*cC**8*sA**8 - 80*cA**8*cB**14*cC**8*sA**8*sB**2 - 344*cA**8*cB**12*cC**8*sA**8*sB**4 + 464*cA**8*cB**10*cC**8*sA**8*sB**6 + 1316*cA**8*cB**8*cC**8*sA**8*sB**8 + 464*cA**8*cB**6*cC**8*sA**8*sB**10 - 344*cA**8*cB**4*cC**8*sA**8*sB**12 - 80*cA**8*cB**2*cC**8*sA**8*sB**14 + 70*cA**8*cC**8*sA**8*sB**16 + 56*cA**6*cB**16*cC**10*sA**10 - 32*cA**6*cB**14*cC**10*sA**10*sB**2 - 288*cA**6*cB**12*cC**10*sA**10*sB**4 + 32*cA**6*cB**10*cC**10*sA**10*sB**6 + 464*cA**6*cB**8*cC**10*sA**10*sB**8 + 32*cA**6*cB**6*cC**10*sA**10*sB**10 - 288*cA**6*cB**4*cC**10*sA**10*sB**12 - 32*cA**6*cB**2*cC**10*sA**10*sB**14 + 56*cA**6*cC**10*sA**10*sB**16 + 28*cA**4*cB**16*cC**12*sA**12 + 32*cA**4*cB**14*cC**12*sA**12*sB**2 - 112*cA**4*cB**12*cC**12*sA**12*sB**4 - 288*cA**4*cB**10*cC**12*sA**12*sB**6 - 344*cA**4*cB**8*cC**12*sA**12*sB**8 - 288*cA**4*cB**6*cC**12*sA**12*sB**10 - 112*cA**4*cB**4*cC**12*sA**12*sB**12 + 32*cA**4*cB**2*cC**12*sA**12*sB**14 + 28*cA**4*cC**12*sA**12*sB**16 + 8*cA**2*cB**16*cC**14*sA**14 + 32*cA**2*cB**14*cC**14*sA**14*sB**2 + 32*cA**2*cB**12*cC**14*sA**14*sB**4 - 32*cA**2*cB**10*cC**14*sA**14*sB**6 - 80*cA**2*cB**8*cC**14*sA**14*sB**8 - 32*cA**2*cB**6*cC**14*sA**14*sB**10 + 32*cA**2*cB**4*cC**14*sA**14*sB**12 + 32*cA**2*cB**2*cC**14*sA**14*sB**14 + 8*cA**2*cC**14*sA**14*sB**16 + cB**16*cC**16*sA**16 + 8*cB**14*cC**16*sA**16*sB**2 + 28*cB**12*cC**16*sA**16*sB**4 + 56*cB**10*cC**16*sA**16*sB**6 + 70*cB**8*cC**16*sA**16*sB**8 + 56*cB**6*cC**16*sA**16*sB**10 + 28*cB**4*cC**16*sA**16*sB**12 + 8*cB**2*cC**16*sA**16*sB**14 + cC**16*sA**16*sB**16'

In general there are limitations on giving radical expressions for the roots of polynomials of degree 5 or more except in certain special cases due to the Abel-Ruffini theorem. If the coefficients are numeric rather than symbolic then we can work around those with things like RootOf. In this case though I doubt that it is possible to express the solutions of this system of equations using anything that you might recognise as an explicit closed form.

If you substitute numeric values for the symbolic unknowns then it should be possible to get approximate roots numerically (e.g. with nroots) or it might become possible to get exact solutions:

In [62]: p_numeric = p.subs({cA:0,sA:1,cB:1/sqrt(2),sB:1/sqrt(2),cC:1})

In [63]: p_numeric
Out[63]: 
                                       8       6       4     2      
  16       14       12       10   35⋅sE    7⋅sE    7⋅sE    sE     1 
sE   - 4⋅sE   + 7⋅sE   - 7⋅sE   + ────── - ───── + ───── - ─── + ───
                                    8        4       16     16   256

In [64]: roots(p_numeric)
Out[64]: 
⎧-√2      √2   ⎫
⎨────: 8, ──: 8⎬
⎩ 2       2    ⎭

In [69]: nroots(p_numeric, n=8, maxsteps=500)
Out[69]: 
[-0.70710678, 0.70710677, 0.70710679, -0.70710679 + 2.4757313e-9⋅ⅈ, -0.70710679 - 2.9650933e-9⋅ⅈ, -
0.70710678 + 4.9626324e-9⋅ⅈ, -0.70710678 - 6.2970508e-9⋅ⅈ, -0.70710678 + 5.3076725e-9⋅ⅈ, -0.7071067
8 - 5.388738e-9⋅ⅈ, -0.70710678 + 2.7281099e-9⋅ⅈ, 0.70710678 + 5.5825484e-9⋅ⅈ, 0.70710678 - 6.328749
4e-9⋅ⅈ, 0.70710678 + 7.8957797e-9⋅ⅈ, 0.70710678 - 8.1564975e-9⋅ⅈ, 0.70710679 + 5.6265597e-9⋅ⅈ, 0.70
710679 - 5.3161607e-9⋅ⅈ]

Note here that nsolve has been unable to determine that there are repeated roots or whether or not the roots are real or have a small imaginary part. This is because finding the roots of a polynomial can be very ill-conditioned especially where there are many repeated roots (possibly the values that I chose to substitute are somewhat degenerate).

Most likely though if you are going to substitute values then it is better to do that at the start and e.g. just call solve/nsolve with the substituted values.

https://en.wikipedia.org/wiki/Resultant

https://en.wikipedia.org/wiki/Abel%E2%80%93Ruffini_theorem

Answer 2

I always like to rewrite trig functions in terms of exponential functions and solve for them, e.g.

>>> from sympy.abc import a
>>> sin(a).rewrite(exp).subs(exp(I*a), x)
-I*(x - 1/x)/2
>>> solve(_)
[-1, 1]

So when exp(I*a) = 1 or -1 you have a solution; that occurs when a is a multiple of pi.

So let's do that with your equations.

If solutions must be real then set up symbols that are real

var('x y',real=True)
con = numbered_symbols('c',real=True)
reps = {} # this will hold replacements that were made

rewrite as expressions that equal 0

eqs = Tuple(*[i.rewrite(Add) for i in (eq1,eq2)])
econ = eqs.replace(
    lambda x: x.is_Function and not x.has(D,E),
    lambda x: reps.setdefault(x, next(con)))
xy = [
    i.rewrite(exp  #rewrite trig in terms of exp
    ).subs(exp(I*D), x).subs(exp(I*E), y  # exp(I*D) -> x, exp(I*E)->y
    ).as_numer_denom()[0].expand() # expand and keep only the numerator
    for i in econ]

split each equation into real and imaginary parts

r,i = zip(*[i.as_real_imag() for i in xy])

If the imaginary parts must be zero then (serendipitously?) those two parts are easy to solve:

solve(i,x,y)
[(-1, y), (0, y), (1, y), (x, 0)]

So when y = exp(I*E) = 0 or +/-1 or x = exp(I*D) = 0 the imaginary parts are 0. Since there are no terms in only x or y we recognize x = 0 or y = 0 as a trivial solution. The case of y = 1 will only make the real parts both zero for special values of the constants as you can check by substituting that value of y into r.

Now if you don't want to consider only real solutions, then the equations will be simultaneously true when the real and imaginary parts are equal:

>>> rdif = factor(r[0]-r[1]); rdif
2*x*y*(2*c0*c5*x**3*y**2 - 2*c0*c5*x**3 - 2*c0*c5*x*y**2 + 2*c0*c5*x +
c1*x**2*y**2 - c1*x**2 - c1*y**2 + c1 + 4*c2*x*y - 4*c3*x**3*y - 4*c3*x*y)
>>> idif = factor(i[0]-i[1]); idif
2*x*y*(-c0 + 2*c1*c5*x)*(x - 1)*(x + 1)*(y**2 + 1)

The imaginary parts are easy to solve for x and/or y and those values can be substituted into rdif and the other value can be solved for. Or SymPy can solve the two equations for you and the solution is not outrageous: Tuple(*solve(idif, rdif, x, y)).count_ops() == 1353. The values for x and y define the corresponding values of exp(I*D) and exp(I*E).

I am not sure what is lost by assuming x and y are real...but hopefully you can see some ways that SymPy can help you restate a problem in a more tractable way.