# How to solve this trigonometric system of equations numerically?

How can the following trigonometric system of equations be solved numerically?

sys= {
1 - Cos[theta1] - 2/3 Sin[phi + π/6] == 0.227746,
h + 2/3 Cos[phi + π/6] - Sin[theta1] == -0.714585,
-1. - Cos[theta2] - 2/3 Sin[phi - π/6] == (3 Cos[psi2])/4,
2/5 + h + 2/3 Cos[phi - π/6] - Sin[theta2] == (3 Sin[psi2])/4,
0.635187 Cos[phi - π/6 - psi2] + 2/3 Cos[1.78586 + phi] Sin[psi2] == 0
}

NSolve does not return a solution even after waiting for a long time:

NSolve[N[sys], {phi, h, theta1, theta2, psi2}]
-
Have a look at FindRoot for non algebraic equations. – b.gatessucks Dec 14 '12 at 10:09
Thank you, this is instantaneous. Although there may be solutions with a very little convergence area, I will use this with a grid-like initialization. Make this comment an answer, I will mark it – Fabio Dalla Libera Dec 14 '12 at 11:29

You can get solutions from NSolve by algebraicizing the system. It works roughly as follows. Expand all trigs, add appropriate trig identities, rename the variables, and the trigonometry disappears for a while. Note that this replacement rule relies on there being no variables that begin with either c or s (more correctly, if there are we will not be able to reverse the replacements as straightforwardly as we do it below).

ruls = {Cos[t_] :> ToExpression[StringJoin["c", ToString[t]]],
Sin[t_] :> ToExpression[StringJoin["s", ToString[t]]]};

sys = TrigExpand[{1 - Cos[theta1] - 2/3 Sin[phi + \[Pi]/6] - 0.227746,
h + 2/3 Cos[phi + \[Pi]/6] - Sin[theta1] - 0.714585, -1. -
Cos[theta2] - 2/3 Sin[phi - \[Pi]/6] - (3 Cos[psi2])/4,
2/5 + h + 2/3 Cos[phi - \[Pi]/6] - Sin[theta2] - (3 Sin[psi2])/4,
0.635187 Cos[phi - \[Pi]/6 - psi2] +
2/3 Cos[1.78586 + phi] Sin[psi2] - 0}];

trigvars = Union[Cases[Variables[sys], (Sin | Cos)[xx_] :> xx]];
trigrels = Map[Sin[#]^2 + Cos[#]^2 - 1 &, trigvars];
fullsys = Join[sys, trigrels] /. ruls

Here is what the system has become.

(* {0.772254 - cphi/3 - ctheta1 - sphi/Sqrt[3], -0.714585 +
cphi/Sqrt[3] + h - sphi/3 - stheta1, -1. + cphi/3 - (3 cpsi2)/4 -
ctheta2 - sphi/Sqrt[3],
2/5 + cphi/Sqrt[3] + h + sphi/3 - (3 spsi2)/4 - stheta2,
0.550088 cphi cpsi2 + 0.317594 cpsi2 sphi - 0.459867 cphi spsi2 -
0.10122 sphi spsi2, -1 + cphi^2 + sphi^2, -1 + cpsi2^2 +
spsi2^2, -1 + ctheta1^2 + stheta1^2, -1 + ctheta2^2 + stheta2^2} *)

We solve it.

sols = NSolve[fullsys];

Length[sols]

(* 20 *)

I would think only real solutions are of interest.

realsols = Select[sols, FreeQ[#, Complex] &]

(* {{h -> -0.418399, ctheta1 -> 0.879946, ctheta2 -> -0.970738,
stheta1 -> -0.475073, stheta2 -> -0.240142, cpsi2 -> 0.792469,
spsi2 -> 0.609913, cphi -> 0.773882,
sphi -> -0.63333}, {h -> -0.132988, ctheta1 -> 0.977023,
ctheta2 -> -0.999392, stheta1 -> -0.213133, stheta2 -> 0.034855,
cpsi2 -> 0.868291, spsi2 -> 0.496056, cphi -> 0.670585,
sphi -> -0.741833}, {h -> 0.0905231, ctheta1 -> 0.999997,
ctheta2 -> 0.325521, stheta1 -> 0.00249828, stheta2 -> 0.945535,
cpsi2 -> -0.892043, spsi2 -> -0.45195, cphi -> 0.643118,
sphi -> -0.765767}, {h -> -0.710519, ctheta1 -> 0.300891,
ctheta2 -> -0.371059, stheta1 -> -0.953659, stheta2 -> 0.928609,
cpsi2 -> -0.608451, spsi2 -> -0.793591, cphi -> 0.965949,
sphi -> 0.258734}} *)

Reversing these trigs to solve for the angles might be done as follows. I remove sines because otherwise I get systems that are overdetermined and, due to use of approximate numbers, (slightly) inconsistent.

trigs = Cases[Variables[sys], (Sin | Cos)[xx_]];
actualvars = Join[trigvars, Complement[Variables[sys], trigs]];

(* {phi, psi2, theta1, theta2, h} *)

varsnosines = Select[Variables[sys], FreeQ[#, Sin] &]

(* {h, Cos[phi], Cos[psi2], Cos[theta1], Cos[theta2]} *)

Here is a way to reverse the trig replacements so we can take our real solutions and replace the altered variables with the original trigs.

revrules = Map[Module[{str = ToString[#]},
If[StringMatchQ[str, "c" ~~ ___], # ->
Cos[ToExpression[StringDrop[str, 1]]],
If[
StringMatchQ[str, "s" ~~ ___], # ->
Sin[ToExpression[StringDrop[str, 1]]]]
]] &, Variables[fullsys]] /. Null :> Sequence[];

(* {cphi -> Cos[phi], cpsi2 -> Cos[psi2], ctheta1 -> Cos[theta1],
ctheta2 -> Cos[theta2], sphi -> Sin[phi], spsi2 -> Sin[psi2],
stheta1 -> Sin[theta1], stheta2 -> Sin[theta2]} *)

We use this to create new systems that are, in a sense, presolved.

newsystems =
Map[Thread[varsnosines - (varsnosines /. #)] &, realsols /. revrules]

(* {{0.418399 + h, -0.773882 + Cos[phi], -0.792469 +
Cos[psi2], -0.879946 + Cos[theta1],
0.970738 + Cos[theta2]}, {0.132988 + h, -0.670585 +
Cos[phi], -0.868291 + Cos[psi2], -0.977023 + Cos[theta1],
0.999392 + Cos[theta2]}, {-0.0905231 + h, -0.643118 + Cos[phi],
0.892043 + Cos[psi2], -0.999997 + Cos[theta1], -0.325521 +
Cos[theta2]}, {0.710519 + h, -0.965949 + Cos[phi],
0.608451 + Cos[psi2], -0.300891 + Cos[theta1],
0.371059 + Cos[theta2]}} *)

Let's solve the first of these for the actual variables of interest.

In[390]:= Solve[newsystems[[1]] == 0, actualvars]

During evaluation of In[390]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. >>

Out[390]= {{phi -> -0.685848, psi2 -> -0.655951, theta1 -> -0.495047, theta2 -> -2.89908, h -> -0.418399}, {phi -> -0.685848, psi2 -> -0.655951, theta1 -> -0.495047, theta2 -> 2.89908, h -> -0.418399}, {phi -> -0.685848, psi2 -> -0.655951, theta1 -> 0.495047, theta2 -> -2.89908, h -> -0.418399}, {phi -> -0.685848, psi2 -> -0.655951, theta1 -> 0.495047, theta2 -> 2.89908, h -> -0.418399}, {phi -> -0.685848, psi2 -> 0.655951, theta1 -> -0.495047, theta2 -> -2.89908, h -> -0.418399}, {phi -> -0.685848, psi2 -> 0.655951, theta1 -> -0.495047, theta2 -> 2.89908, h -> -0.418399}, {phi -> -0.685848, psi2 -> 0.655951, theta1 -> 0.495047, theta2 -> -2.89908, h -> -0.418399}, {phi -> -0.685848, psi2 -> 0.655951, theta1 -> 0.495047, theta2 -> 2.89908, h -> -0.418399}, {phi -> 0.685848, psi2 -> -0.655951, theta1 -> -0.495047, theta2 -> -2.89908, h -> -0.418399}, {phi -> 0.685848, psi2 -> -0.655951, theta1 -> -0.495047, theta2 -> 2.89908, h -> -0.418399}, {phi -> 0.685848, psi2 -> -0.655951, theta1 -> 0.495047, theta2 -> -2.89908, h -> -0.418399}, {phi -> 0.685848, psi2 -> -0.655951, theta1 -> 0.495047, theta2 -> 2.89908, h -> -0.418399}, {phi -> 0.685848, psi2 -> 0.655951, theta1 -> -0.495047, theta2 -> -2.89908, h -> -0.418399}, {phi -> 0.685848, psi2 -> 0.655951, theta1 -> -0.495047, theta2 -> 2.89908, h -> -0.418399}, {phi -> 0.685848, psi2 -> 0.655951, theta1 -> 0.495047, theta2 -> -2.89908, h -> -0.418399}, {phi -> 0.685848, psi2 -> 0.655951, theta1 -> 0.495047, theta2 -> 2.89908, h -> -0.418399}}

Okay, they are pretty much all the same up to sign. Maybe take the last one since it is mostly positive (h is always negative; only the trig vars vary).

We can do the same with all of the new systems. I'll extract the last solution from each as being a representative.

Quiet[
usefulsols = Map[Last[Solve[# == 0, actualvars]] &, newsystems]]

(* {{phi -> 0.685848, psi2 -> 0.655951, theta1 -> 0.495047,
theta2 -> 2.89908, h -> -0.418399}, {phi -> 0.8358, psi2 -> 0.51905,
theta1 -> 0.21478, theta2 -> 3.10673,
h -> -0.132988}, {phi -> 0.872233, psi2 -> 2.67264,
theta1 -> 0.00249829, theta2 -> 1.23923,
h -> 0.0905231}, {phi -> 0.261712, psi2 -> 2.2249,
theta1 -> 1.26517, theta2 -> 1.95095, h -> -0.710519}} *)
-
I am curious, does NSolve attempt the same (very straightforward) method that you described here? – Szabolcs Dec 14 '12 at 16:55
@Szabolcs To paraphrase Gertude Stein, there is no here here (I think you missed a link). – Daniel Lichtblau Dec 14 '12 at 17:19
I meant the method you describe in this answer above: to convert the system of trig equations to a system of polynomial equations, and solve that. It seems like something that could be automated, so I was wondering if this is what NSolve tries – Szabolcs Dec 14 '12 at 21:39
@Szabolcs No, NSolve does not do any of that. – Daniel Lichtblau Dec 14 '12 at 21:48
This is a great answer, thank you – Fabio Dalla Libera Dec 15 '12 at 2:01