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