Re: Help with non-linear system of equations

*To*: mathgroup at smc.vnet.net*Subject*: [mg115816] Re: Help with non-linear system of equations*From*: Manjinder Benning <manjbenning at gmail.com>*Date*: Fri, 21 Jan 2011 04:34:01 -0500 (EST)

thanks Daniel i wanted to confirm that i am normalizing my ellipse correctly Take, x^2/a^2 + y^2/b^2 = 1 Where the a and b are half my major and minor axes lengths. Then i define a1 = b^2/a^2 b1 = b^2 this gives me the normalized form that you used: ellipse = a1*x^2 + y^2 - b1 Is this how you meant? I ran some tests using your proposed solution and found that it doesnt seem to work. Ex: I collected~20 points and fit an ellipse to it. Major Axis Length: 137.6 Minor Axis Length: 123 Using. Solve[(polys /. {b -> 3782, a -> 0.8}) == 0] // N Out[31]{{cossqr->-0.01838,xc->-30.9734,zp->0.-7.38965 I,xp->-30.5262},{cossqr->-0.01838,xc->-30.9734,zp->0.+7.38965 I,xp->-30.5262},{cossqr->-0.01838,xc->30.9734,zp->0.-7.38965 I,xp->30.5262},{cossqr->-0.01838,xc->30.9734,zp->0.+7.38965 I,xp->30.5262},{cossqr->0.0179834,xc->-30.5262,zp->-7.44358,xp->-30.9734},{cossqr->0.0179834,xc->-30.5262,zp->7.44358,xp->-30.9734},{cossqr->0.0179834,xc->30.5262,zp->-7.44358,xp->30.9734},{cossqr->0.0179834,xc->30.5262,zp->7.44358,xp->30.9734}} On Tue, Jan 18, 2011 at 2:05 PM, Daniel Lichtblau <danl at wolfram.com> wrote: > 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 > -- www.thepowerofyang.com

**Re: Do I need MathLink to run finite-difference fast enough for Manipulate?**

**Re: Help with non-linear system of equations**

**Re: Compilation options question**

**Re: Help with non-linear system of equations**