Re: Help with non-linear system of equations (Trigonometric and Polynomial)
- To: mathgroup at smc.vnet.net
- Subject: [mg115730] Re: Help with non-linear system of equations (Trigonometric and Polynomial)
- From: Daniel Lichtblau <danl at wolfram.com>
- 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
that.
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 =
Numerator[
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
example.
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