How to solve this trigonometric system of equations numerically?

Question

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}]

Answer 1

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

Answer 2

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
     }
   }
 ]