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