Re: Help with non-linear system of equations

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

```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

```

• Prev by Date: Re: Help with non-linear system of equations
• Next by Date: Re: Help with non-linear system of equations
• Previous by thread: Re: Help with non-linear system of equations
• Next by thread: Re: Help with non-linear system of equations