Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive

MathGroup Archive 2011

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Help with non-linear system of equations (Trigonometric and Polynomial)

  • To: mathgroup at
  • Subject: [mg115730] Re: Help with non-linear system of equations (Trigonometric and Polynomial)
  • From: Daniel Lichtblau <danl at>
  • Date: Wed, 19 Jan 2011 05:28:20 -0500 (EST)

Manjinder Benning wrote:
> Hi i am trying to do some 3D modeling by finding a 6 unknown 6 equation
> simultaneous system. I am using the equation:
> Cos(a) = A.B/|A||B|
> where A and B are lines (vectors) in 3 space.
> a is the angle between the lines.
> When i run the calculation, Mathematica just runs forever and doesn't come
> to a solution. It would be great if a closed form solution could be found.
> I am understanding that it is difficult to solve systems with polynomial and
> trig parts?
> Description of Problem:
> So i have 3 points in 3-space.
> Pl (Xl, Yl, Zl) - A constant point in the air
> Pc(Xc, Yc, Zc) - A constant point on the ground (Z=0)
> Pe's (Xe_n, Ye_n, Ze_n) - A series of points that lay on the edges of an
> ellipse co-planar to Z=0 (which are known experimentally)
> A - The angle between the line (Pl -> Pc) and each (Pl -> Pe_n)
> Due to the geometry of the problem (a ellipse) angle A is constant.
> Are there any methods i can use to make it easier to solve (ie bounding the
> problem) or am i dreaming that i will find a closed form solution?
> ...any help would be appreciated:
> Heres the code:
> eq1 =Cos[A] ==  ((Xc - Xl)*(Xe1-Xl) + (Yc-Yl)*(Ye1-Yl) + (-Zl)*(-Zl))/Sqrt[(
> (Xc - Xl)^2 + (Yc-Yl)^2+(-Zl)^2 ) * ((Xe1-Xl)^2 +(Ye1-Yl)^2+(-Zl)^2)]
> eq2 =Cos[A] ==  ((Xc - Xl)*(Xe2-Xl) + (Yc-Yl)*(Ye2-Yl) + (-Zl)*(-Zl))/Sqrt[(
> (Xc - Xl)^2 + (Yc-Yl)^2+(-Zl)^2 ) * ((Xe2-Xl)^2 +(Ye2-Yl)^2+(-Zl)^2)]
> eq3 =Cos[A] ==  ((Xc - Xl)*(Xe3-Xl) + (Yc-Yl)*(Ye3-Yl) + (-Zl)*(-Zl))/Sqrt[(
> (Xc - Xl)^2 + (Yc-Yl)^2+(-Zl)^2 ) * ((Xe3-Xl)^2 +(Ye3-Yl)^2+(-Zl)^2)]
> eq4 =Cos[A] ==  ((Xc - Xl)*(Xe4-Xl) + (Yc-Yl)*(Ye4-Yl) + (-Zl)*(-Zl))/Sqrt[(
> (Xc - Xl)^2 + (Yc-Yl)^2+(-Zl)^2 ) * ((Xe4-Xl)^2 +(Ye4-Yl)^2+(-Zl)^2)]
> eq5 =Cos[A] ==  ((Xc - Xl)*(Xe5-Xl) + (Yc-Yl)*(Ye5-Yl) + (-Zl)*(-Zl))/Sqrt[(
> (Xc - Xl)^2 + (Yc-Yl)^2+(-Zl)^2 ) * ((Xe5-Xl)^2 +(Ye5-Yl)^2+(-Zl)^2)]
> eq6 =Cos[A] ==  ((Xc - Xl)*(Xe6-Xl) + (Yc-Yl)*(Ye6-Yl) + (-Zl)*(-Zl))/Sqrt[(
> (Xc - Xl)^2 + (Yc-Yl)^2+(-Zl)^2 ) * ((Xe6-Xl)^2 +(Ye6-Yl)^2+(-Zl)^2)]
> Solve[{eq1,eq2,eq3,eq4,eq5,
> eq6},{A,Xl,Yl,Zl,Xc,Yc}]
> much gratitude to all things Mathematica!
> Manjinder.

Sjoerd de Vries gave a good answer but there are some fine points I had 
wanted to address. Please pardon a certain amount of replication of what 
he mentioned.

First is that your system is not quite correct. It is overdetermined in 
the sense that you give six points on an ellipse, and underdetermined in 
that it will not be able to solve for your six parameters. The reason is 
of course that once you have the point coordinates, your angle is forced 
by e.g. the Law of Cosines. I assert this boldly so nobody will realize 
it took me two days to locate my missing equation. Anyway, you now can 
form a sixth equation using the cosine law.

How do we know the original formulation is really overdetermined? 
because the geometry of constant angle from a line does in fact 
determine a quadratic curve (a conic section in the plane), and five 
points suffice to uniquely identify such a conic.

Next item is that in some sense a closed form solution must exist. But 
it is likely to be horrendously big and not numerically well behaved 
when you plug in actual values for your points. Since your data is 
experimental there seems to be little point in not opting for a more 
numeric approach. Moreover such an approach might behave better in the 
situation of having an overdetermined ellipse (as de Vries notes, you 
can just obtain said ellipse via a best fit to the six points).

It simplifies matters considerably to move the ellipse so that the major 
axis is the x axis. You can also center it at the origin, though that 
does not really affect the equations we will use. The equations needed 
for translating and rotating an ellipse are easy to set up and I'll omit 

We really want to work with polynomials, so we'll take some squares as 
needed. Also notice that the angle only appears as as cosine, and we'll 
end up solving for cos-squared. So there really is no transcendental 
function until we're done with the algebra, at which point an arccosine 
will be needed to recover the angle.

We'll assume we normalize the ellipse so that it has the form 
a*x^2+y^2-b==0, with 0<a<=1. Some geometric consideration indicates that 
the segment from ground point to air point must project onto the major 
axis, so we do not have yc or yp to bother with. So we have parameters 
xc, xp, zp, and cossqr.

Here is a symbolic formulation.

groundpt = {xc, yc, 0};
airpt = {xp, yp, zp};
groundpt = {xc, 0, 0};
airpt = {xp, 0, zp};
ellipse = a*x^2 + y^2 - b;
cossquared[v1_, v2_] := (v1.v2)^2/(v1.v1*v2.v2)

Here we get an expression to the effect that the cosine squared of angle 
between segment from air point to ellipse point and segment from air 
point to ground point, is some constant value. We have, by set-up, 
assumed the major axis is the x axis.

In[132]:= expr =
   Together[cossqr - cossquared[{x, y, 0} - airpt, groundpt - airpt]]]

Out[132]= -x^2 xc^2 + cossqr x^2 xc^2 + 2 x^2 xc xp -
  2 cossqr x^2 xc xp + 2 x xc^2 xp - 2 cossqr x xc^2 xp - x^2 xp^2 +
  cossqr x^2 xp^2 - 4 x xc xp^2 + 4 cossqr x xc xp^2 - xc^2 xp^2 +
  cossqr xc^2 xp^2 + 2 x xp^3 - 2 cossqr x xp^3 + 2 xc xp^3 -
  2 cossqr xc xp^3 - xp^4 + cossqr xp^4 + cossqr xc^2 y^2 -
  2 cossqr xc xp y^2 + cossqr xp^2 y^2 + cossqr x^2 zp^2 -
  2 x xc zp^2 + cossqr xc^2 zp^2 + 2 x xp zp^2 - 2 cossqr x xp zp^2 +
  2 xc xp zp^2 - 2 cossqr xc xp zp^2 - 2 xp^2 zp^2 +
  2 cossqr xp^2 zp^2 + cossqr y^2 zp^2 - zp^4 + cossqr zp^4

Notice the above is a quadratic in x and y, that is, the locus in 
question does in fact lie on a conic.

To get equations for the point and angle parameters, in terms of our 
ellipse parameters, we equate terms in the ellipse.

In[155]:= polys =
  Flatten[CoefficientList[expr, {x, y}]] -
    Flatten[CoefficientList[ellipse, {x, y}]] /. 0 :> Sequence[]

Out[155]= {b - xc^2 xp^2 + cossqr xc^2 xp^2 + 2 xc xp^3 -
   2 cossqr xc xp^3 - xp^4 + cossqr xp^4 + cossqr xc^2 zp^2 +
   2 xc xp zp^2 - 2 cossqr xc xp zp^2 - 2 xp^2 zp^2 +
   2 cossqr xp^2 zp^2 - zp^4 + cossqr zp^4, -1 + cossqr xc^2 -
   2 cossqr xc xp + cossqr xp^2 + cossqr zp^2,
  2 xc^2 xp - 2 cossqr xc^2 xp - 4 xc xp^2 + 4 cossqr xc xp^2 +
   2 xp^3 - 2 cossqr xp^3 - 2 xc zp^2 + 2 xp zp^2 -
   2 cossqr xp zp^2, -a - xc^2 + cossqr xc^2 + 2 xc xp -
   2 cossqr xc xp - xp^2 + cossqr xp^2 + cossqr zp^2}

Say we give numeric values for the ellipse parameters. Then it is a 
simple matter to find the point and angle parameters. Here is a simple 

In[156]:= Solve[(polys /. {b -> 1, a -> 1/9}) == 0] // N

Out[156]= {{cossqr -> 0.779803, xc -> 2.39604, zp -> -0.627285,
   xp -> 3.33885}, {cossqr -> 0.779803, xc -> 2.39604, zp -> 0.627285,
   xp -> 3.33885}, {cossqr -> 0.779803, xc -> -2.39604,
   zp -> -0.627285, xp -> -3.33885}, {cossqr -> 0.779803,
   xc -> -2.39604, zp -> 0.627285, xp -> -3.33885}, {cossqr -> 1.64877,
    xc -> -3.33885, zp -> 0. - 0.53139 I,
   xp -> -2.39604}, {cossqr -> 1.64877, xc -> -3.33885,
   zp -> 0.+ 0.53139 I, xp -> -2.39604},
  {cossqr -> 1.64877, xc -> 3.33885,
   zp -> 0. - 0.53139 I,
   xp -> 2.39604}, {cossqr -> 1.64877, xc -> 3.33885,
   zp -> 0. + 0.53139 I, xp -> 2.39604}}

Complex solutions (where the cosine squared may be other than between 0 
and 1) get discarded. You then can decide what you want to use. 
Presumably you want the z coordinate of the air point to be positive.

Another subtlety here is that we have the square of the cosine. Should 
we take the cosine to be positive or negative? This depends on the 
placement of xp. Geometric considerations show that xc will be inside 
the ellipse. If xp is such that the air point projects (down) to be also 
inside the ellipse then the angle is obtuse, so take the negative for 
the cosine. Else it is acute and so the cosine is positive.

By the way, while the above use of Solve (or NSolve) seems easy enough, 
a different formulation with FindRoot seemed unable to get solutions 
even when starting near one. I do not know if this was bad setup on my 
part or numeric instability. In any case, what this indicates to me is 
that the steps of finding the ellipse from six points, 
translating/rotating to put into the above canonical form, and solving 
for {xc,xp,zp,cossqr} are the way to do this. In contrast, setting up 
five equations based on five ellipse points, augmenting with an equation 
derived from the Law of Cosines, and solving from there seems not to 
work: FindRoot gives up and NSolve hangs.

For the morbidly curious, the extra eqwuation based on LoC is 
constructed as below.

newpoly = With[
   {v1 = groundpt - pts[[1]], v2 = groundpt - airpt,
    v3 = pts[[1]] - airpt},
   (v1.v1 - v2.v2 - v3.v3)^2 - 4*v2.v2*v3.v3*cossqr]

Daniel Lichtblau
Wolfram Research

  • Prev by Date: Re: Array of arrays of various sizes and compile
  • Next by Date: Re: CUDADot bug?
  • Previous by thread: Re: Help with non-linear system of equations (Trigonometric and Polynomial)
  • Next by thread: How to force numeric output