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