Re: Help with non-linear system of equations
- To: mathgroup at smc.vnet.net
- Subject: [mg115830] Re: Help with non-linear system of equations
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Sat, 22 Jan 2011 03:21:03 -0500 (EST)
Manjinder Benning wrote:
> I guess as the ellipse ----> circle there will be infinite combinations of
> zc and the angle..hmm......
Thank you for answering your own question so I don't have to. I'll
provide a bit of detail.
(1) You are normalizing your ellipse in the way I had in mind.
(2) The circle case is indeed the one where this won't work. Unlike the
others, it has infinitely many solutions. This can only happen on a
measure zero set in your parameter space (the space of all possible
quadratic curves, in this case). Not a big deal, because you can
special-case to recognize when your points lie on a circle.
(3) You really cannot restrict the angle parameter. It is only when you
have a bona fide circle that any angle will work. In all other cases
there are only four (real valued) solutions, and they all involve the
same angle.
(4) What happens as the ellipse gets rounder (approaches a circle) is
that the number of solutions remains the same. The ONLY bad case, in the
sense of number of solutions, is the circle itself, not nearby. But you
are quite correct that you will encounter trouble as you get near to a
circle. It manifests as numerical ill conditioning of the equations of
interest. Offhand I do not have a good solution for that; my suggestion
is to use NSolve with high precision e.g. WorkingPrecision->200. If you
do that, you may need to supply your approximated ellipse with high
precision numbers. Not a great scenario, since the thing is only crudely
approximated via best fit to begin with. But I don't have a better one
in mind at the moment.
Daniel Lichtblau
Wolfram Research
> 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
>>
>
>
>