# 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

Using Manipulate and FindRoot makes it easier to look for roots. I do not know why NSolve can't do it. Hard problem.

Normally one starts with NSolve, then FindRoot. I looked to see what other options are there, and could not find one (since these are not polynomials, else I would tried LinearSolve on them also, but can't since these are trig equations).

Using this below at least, you can look for roots by changing the parameters more easily.

btw, You will get warnings such as these using FindRoot as is on this problem:

FindRoot::jsing: Encountered a singular Jacobian at the point...Try perturbing the initial point(s). >>

And

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances. >>

You can look into ?FindRoot for more information on how to change some of its options to try to avoid these warning if possible.

Manipulate[
TableForm@
FindRoot[sys, {{phi,phi0}, {h, h0},{theta1,theta10},{theta2,theta20},{psi2, psi20}}],

{{phi0, 1, "around which phi?"}, -5, 5, .01,Appearance->"Labeled",ImageSize -> Small},
{{h0, 1, "around which h?"}, -5, 5, .01, Appearance -> "Labeled",ImageSize -> Small},
{{theta10,1,"around which theta1?"},-5, 5,.01,Appearance->"Labeled",ImageSize->Small},
{{theta20,1,"around which theta2?"},-5, 5,.01,Appearance->"Labeled",ImageSize->Small},
{{psi20,1,"around which psi2?"},-5,5,.01,Appearance -> "Labeled", ImageSize -> Small},
ContentSize -> {300, 100},
Initialization :> {sys = {
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
}
}
]


-
Good answer and a nice tip for how to build a tool for exploration in these cases. Two comments: (1) You indicate that you can execute ?NSolve to understand the warning messages - I assume you mean ?FindRoot. (2) NSolve employs hybrid numeric/symbolic schemes that don't work too well with transcendental equations - in fact, the documentation indicates that it "deals primarily with linear and polynomial equations"; that's why it doesn't work. – Mark McClure Dec 14 '12 at 12:55
Yes, I meant to write ?FindRoot, will correct. thanks. – Nasser Dec 14 '12 at 13:31