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

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

score 11 · Accepted Answer · edited Dec 14 '12 at 17:14

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

How to solve this trigonometric system of equations numerically?

1 Answers1

Linked

Related