Re: Help with non-linear system of equations

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

I guess as the ellipse ----> circle there will be infinite combinations of zc and the angle..hmm...... On Thu, Jan 20, 2011 at 7:41 PM, Manjinder Benning <manjbenning at gmail.com>wrote: > opps...accidentally sent that before it was complete:------- > > > 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: > Pp is directly above Pc (only differ by height yp = yc, xp = yc) > I collected ~20 points and fit an ellipse to them. > > Major Axis Length: 137.6 > Minor Axis Length: 135 > Using. > Solve[(polys /. {b -> 4556, a -> 0.96}) == 0] // N > > {{cossqr->-0.0152402,xc->-13.8784,zp->0.-8.10285 > I,xp->-13.6784},{cossqr->-0.0152402,xc->-13.8784,zp->0.+8.10285 > I,xp->-13.6784},{cossqr->-0.0152402,xc->13.8784,zp->0.-8.10285 > I,xp->13.6784},{cossqr->-0.0152402,xc->13.8784,zp->0.+8.10285 > I,xp->13.6784},{cossqr->0.0150024,xc->-13.6784,zp->-8.16187,xp->-13.8784},{cossqr->0.0150024,xc->-13.6784,zp->8.16187,xp->-13.8784},{cossqr->0.0150024,xc->13.6784,zp->-8.16187,xp->13.8784},{cossqr->0.0150024,xc->13.6784,zp->8.16187,xp->13.8784}} > So out of this i get that > xp = 13.87 > xc = 13.67 > zc = 8.16 > cossqr = 0.015 .......angle = 82.9 degrees. > > The actual values are: > xp = ~0 > xc = ~0 > zc = ~ 70 > Angle = ~45 degrees > > The xp and xc seem to be close (error lies in my ellipse fit), however the > height zc and angle are way off. > I guess this solution sort of makes sense that if zc was so low as 8cm then > the angle would have to be as high as 83 degrees. > So it looks like another solution has been found other than the one that > created the 20 ellipse points. > > Is there a way i can bound the angle to be say < 50 degrees? > > > Manjinder > > > 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 > -- www.thepowerofyang.com

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

**Re: adding a keyboard shortcut for double brackets**

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

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