Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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}]
share|improve this question
    
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
add comment

1 Answer

up vote 8 down vote accepted

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}} *)
share|improve this answer
    
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.