Got Questions? Get Answers.
Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)

Subject: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)

From: GAURAV

Date: 29 Apr, 2010 20:49:06

Message: 1 of 18

Hello Eveyone,

I am trying to get an ellipse passing through the four points I have which are not symmetric.
So, I am using the equation (X-Xc)^2/(a^2) + (Y-Yc)^2/(b^2) =1

I have four points which wil go in X and Y of the equation above.
Let the points be (a1,b1), (a2,b2), (a3,b3) and (a4,b4)..

How do I get a,b,Xc and Yc from this??

Please help me out. Really appreciate..

Thanks,
Gaurav

Subject: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)

From: Roger Stafford

Date: 29 Apr, 2010 21:03:04

Message: 2 of 18

"GAURAV " <gsharda@engineering.uiowa.edu> wrote in message <hrcrc2$3me$1@fred.mathworks.com>...
> Hello Eveyone,
>
> I am trying to get an ellipse passing through the four points I have which are not symmetric.
> So, I am using the equation (X-Xc)^2/(a^2) + (Y-Yc)^2/(b^2) =1
>
> I have four points which wil go in X and Y of the equation above.
> Let the points be (a1,b1), (a2,b2), (a3,b3) and (a4,b4)..
>
> How do I get a,b,Xc and Yc from this??
>
> Please help me out. Really appreciate..
>
> Thanks,
> Gaurav

  Are you certain the ellipse you are looking for has its major and minor axes aligned with the x and y axes? If not, a minimum of five points is necessary.

Roger Stafford

Subject: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)

From: GAURAV

Date: 29 Apr, 2010 21:21:04

Message: 3 of 18

Hi Roger,
Thanks for your reply..

Actually its an assumption I am making.
I have two curves(orthogonal to each other) aligned in 3D ( not necessarily symmetrical ).
The height (along Z axis) is the same. So, when I cut a slice for some Z, I have four points (two from each curve).
Now I need to construct a volume from these two views. So from each slice I am trying to fit in an ellipse. I hope this makes sense.

So, now I have four points which pass through ellipse and ellipse should not go out of this bound .
So, I am looking to solve the ellipse equation and try to find out of the a,b,xc and yc required for the ellipse with angle=0.

Can you suggest something?

Regards,
Gaurav.

Subject: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)

From: GAURAV

Date: 29 Apr, 2010 21:22:04

Message: 4 of 18

Hi Roger,
Thanks for your reply..

Actually its an assumption I am making.
I have two curves(orthogonal to each other) aligned in 3D ( not necessarily symmetrical ).
The height (along Z axis) is the same. So, when I cut a slice for some Z, I have four points (two from each curve).
Now I need to construct a volume from these two views. So from each slice I am trying to fit in an ellipse. I hope this makes sense.

So, now I have four points which pass through ellipse and ellipse should not go out of this bound .
So, I am looking to solve the ellipse equation and try to find out of the a,b,xc and yc required for the ellipse with angle=0.

Can you suggest something?

Regards,
Gaurav.

Subject: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)

From: Roger Stafford

Date: 29 Apr, 2010 22:16:04

Message: 5 of 18

"GAURAV " <gsharda@engineering.uiowa.edu> wrote in message <hrct9s$bpa$1@fred.mathworks.com>...
> Hi Roger,
> Thanks for your reply..
>
> Actually its an assumption I am making.
> I have two curves(orthogonal to each other) aligned in 3D ( not necessarily symmetrical ).
> The height (along Z axis) is the same. So, when I cut a slice for some Z, I have four points (two from each curve).
> Now I need to construct a volume from these two views. So from each slice I am trying to fit in an ellipse. I hope this makes sense.
>
> So, now I have four points which pass through ellipse and ellipse should not go out of this bound .
> So, I am looking to solve the ellipse equation and try to find out of the a,b,xc and yc required for the ellipse with angle=0.
>
> Can you suggest something?
>
> Regards,
> Gaurav.

  For the four-point problem you can avoid using one of the non-linear methods by writing your equation in the following linear form. Then use 'svd' to determine the five unknown coefficients:

 A*x^2 + C*y^2 + D*x + E*y + F = 0

Do this:

 x = [a1;a2;a3;a4];
 y = [b1;b2;b3;b4];
 M = [x.^2,y.^2,x,y,ones(4,1)];
 [U,S,V] = svd(M);
 A = V(1,5);
 C = V(2,5);
 D = V(3,5);
 E = V(4,5);
 F = V(5,5);

The fifth column of V has your desired coefficients. These five coefficients should satisfy the above equation.

  By completion of the squares in the equation and dividing by the appropriate value you can then transform it into your original form and thereby determine a, b, xc, and xc.

  By way of explanation, this fifth column is the one corresponding to the zero-valued singular value of the matrix M, which causes it to be a solution to the equation. There is bound to be at least one singular value because the number of rows in M is less that its number of columns, and 'svd' always puts the least singular value in the last column - the fifth in this case.

Roger Stafford

Subject: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)

From: Walter Roberson

Date: 29 Apr, 2010 22:18:27

Message: 6 of 18

GAURAV wrote:
> Hello Eveyone,
>
> I am trying to get an ellipse passing through the four points I have
> which are not symmetric.
> So, I am using the equation (X-Xc)^2/(a^2) + (Y-Yc)^2/(b^2) =1
>
> I have four points which wil go in X and Y of the equation above.
> Let the points be (a1,b1), (a2,b2), (a3,b3) and (a4,b4)..
>
> How do I get a,b,Xc and Yc from this??

Answer 1:

[a = 0, b = 0, Xc = Xc, Yc = Yc]

Answer 2:

[a = (1/4) * 4^(1/2) * ((-a4 * b3 + a4 * b2 - b4 * a2 - a3 * b2 + a3 *
b4 + b3 * a2) * (-a3 * b2 + a4 * b2 + a1 * b3 - a4 * b3 - a1 * b4 + a3 *
b4) * (a3 * b4 - a3 * b2 + a1 * b2 + a2 * b2 - a4 * b2 - a1 * b3 + a4 *
b3 - b4 * a2) * (a3 * b4 - a3 * b2 + a1 * b2 + a2 * b2 - a4 * b2 - b3 *
a2 + a4 * b3 - a1 * b4) / ((-b4 + b3) * (-b4 + b2) * (-b3 + b2) * (-a4 *
a1 * b3 + b2 * a4 * a1 - a1 * a2 * b4 - b2 * a3 * a1 + b4 * a3 * a1 + a1
* a2 * b3 + a3 * a2 * b4 - a4 * a2 * b3 + b2 * a3^2 - b2 * a3 * a2 -
a3^2 * b4 + b3 * a4^2 + b2 * a4 * a2 - b2 * a4^2)))^(1/2),

b = (1/2) * ((-a4 * b3 + a4 * b2 - b4 * a2 - a3 * b2 + a3 * b4 + b3 *
a2) * (-a3 * b2 + a4 * b2 + a1 * b3 - a4 * b3 - a1 * b4 + a3 * b4) * (a3
* b4 - a3 * b2 + a1 * b2 + a2 * b2 - a4 * b2 - a1 * b3 + a4 * b3 - b4 *
a2) * (a3 * b4 - a3 * b2 + a1 * b2 + a2 * b2 - a4 * b2 - b3 * a2 + a4 *
b3 - a1 * b4))^(1/2) / (((-b4 + b3) * a2 + (b4 - b2) * a3 + (-b3 + b2) *
a4) * a1 + ((b4 - b2) * a3 + (-b3 + b2) * a4) * a2 + (-b4 + b2) * a3^2 +
(-b2 + b3) * a4^2),

Xc = (1/2) * a1 + (1/2) * a2,

Yc = (1/2) * (((b3^2 - b4^2) * a2 + (-b2^2 + b4^2) * a3 + (-b3^2 + b2^2)
* a4) * a1 + ((-b2^2 + b4^2) * a3 + (-b3^2 + b2^2) * a4) * a2 + (b2^2 -
b4^2) * a3^2 + (b3^2 - b2^2) * a4^2) / (((-b4 + b3) * a2 + (b4 - b2) *
a3 + (-b3 + b2) * a4) * a1 + ((b4 - b2) * a3 + (-b3 + b2) * a4) * a2 +
(-b4 + b2) * a3^2 + (-b2 + b3) *a4^2)]

Answer 3:

[a = -(1/2) * ((-a4 * b3 + a4 * b2 - b4 * a2 - a3 * b2 + a3 * b4 + b3 *
a2) * (-a3 * b2 + a4 * b2 + a1 * b3 - a4 * b3 - a1 * b4 + a3 * b4) * (a3
* b4 - a3 * b2 + a1 * b2 + a2 * b2 - a4 * b2 - a1 * b3 + a4 * b3 - b4 *
a2) * (a3 * b4 - a3 * b2 + a1 * b2 + a2 * b2 - a4 * b2 - b3 * a2 + a4 *
b3 - a1 * b4) / ((-b4 + b3) * (-b4 + b2) * (-b3 + b2) * (-a4 * a1 * b3 +
b2 * a4 * a1 - a1 * a2 * b4 - b2 * a3 * a1 + b4 * a3 * a1 + a1 * a2 * b3
+ a3 * a2 * b4 - a4 * a2 * b3 + b2 * a3^2 - b2 * a3 * a2 - a3^2 * b4 +
b3 * a4^2 + b2 * a4 * a2 - b2 * a4^2)))^(1/2),

b = (1/2) * ((-a4 * b3 + a4 * b2 - b4 * a2 - a3 * b2 + a3 * b4 + b3 *
a2) * (-a3 * b2 + a4 * b2 + a1 * b3 - a4 * b3 - a1 * b4 + a3 * b4) * (a3
* b4 - a3 * b2 + a1 * b2 + a2 * b2 - a4 * b2 - a1 * b3 + a4 * b3 - b4 *
a2) * (a3 * b4 - a3 * b2 + a1 * b2 + a2 * b2 - a4 * b2 - b3 * a2 + a4 *
b3 - a1 * b4))^(1/2) / (((-b4 + b3) * a2 + (b4 - b2) * a3 + (-b3 + b2) *
a4) * a1 + ((b4 - b2) * a3 + (-b3 + b2) * a4) * a2 + (-b4 + b2) * a3^2 +
(-b2 + b3) * a4^2),

Xc = (1/2) * a1 + (1/2) * a2,

Yc = (1/2) * (((b3^2 - b4^2) * a2 + (-b2^2 +b4^2) * a3 + (-b3^2 + b2^2)
* a4) * a1 + ((-b2^2 + b4^2) * a3 + (-b3^2 + b2^2) * a4) * a2 + (b2^2 -
b4^2) * a3^2 + (b3^2 - b2^2) * a4^2) / (((-b4 + b3) * a2 + (b4 - b2) *
a3 + (-b3 + b2) * a4) * a1 + ((b4 - b2) * a3 + (-b3 + b2) * a4) * a2 +
(-b4 + b2) * a3^2 + (-b2 + b3) * a4^2)]

Answer 4:

[a = (1/4) * 4^(1/2) * ((-a4 * b3 + a4 * b2 - b4 * a2 - a3 * b2 + a3 *
b4 + b3 * a2) * (-a3 * b2 + a4 * b2 + a1 * b3-a4 * b3-a1 * b4 + a3 *
b4) * (a3 * b4 - a3 * b2 + a1 * b2 + a2 * b2 - a4 * b2 - a1 * b3 + a4 *
b3 - b4 * a2) * (a3 * b4 - a3 * b2 + a1 * b2 + a2 * b2 - a4 * b2 - b3 *
a2 + a4 * b3 - a1 * b4)/((-b4 + b3) * (-b4 + b2) * (-b3 + b2) * (-a4 *
a1 * b3 + b2 * a4 * a1 - a1 * a2 * b4 - b2 * a3 * a1 + b4 * a3 * a1 + a1
* a2 * b3 + a3 * a2 * b4 - a4 * a2 * b3 + b2 * a3^2 - b2 * a3 * a2 -
a3^2 * b4 + b3 * a4^2 + b2 * a4 * a2 - b2 * a4^2)))^(1/2),

b = -(1/2) * ((-a4 * b3 + a4 * b2 - b4 * a2 - a3 * b2 + a3 * b4 + b3 *
a2) * (-a3 * b2 + a4 * b2 + a1 * b3 - a4 * b3 - a1 * b4 + a3 * b4) * (a3
* b4 - a3 * b2 + a1 * b2 + a2 * b2 - a4 * b2 - a1 * b3 + a4 * b3 - b4 *
a2) * (a3 * b4 - a3 * b2 + a1 * b2 + a2 * b2 - a4 * b2 - b3 * a2 + a4 *
b3 - a1 * b4))^(1/2)/(((-b4 + b3) * a2 + (b4 - b2) * a3 + (-b3 + b2) *
a4) * a1 + ((b4 - b2) * a3 + (-b3 + b2) * a4) * a2 + (-b4 + b2) * a3^2 +
(-b2 + b3) * a4^2),

Xc = (1/2) * a1 + (1/2) * a2,

Yc = (1/2) * (((b3^2 - b4^2) * a2 + (-b2^2 + b4^2) * a3 + (-b3^2 + b2^2)
* a4) * a1 + ((-b2^2 + b4^2) * a3 + (-b3^2 + b2^2) * a4) * a2 + (b2^2 -
b4^2) * a3^2 + (b3^2 - b2^2) * a4^2)/(((-b4 + b3) * a2 + (b4 - b2) * a3
+ (-b3 + b2) * a4) * a1 + ((b4 - b2) * a3 + (-b3 + b2) * a4) * a2 + (-b4
+ b2) * a3^2 + (-b2 + b3) * a4^2)]

Subject: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)

From: Roger Stafford

Date: 29 Apr, 2010 22:41:27

Message: 7 of 18

Walter Roberson <roberson@hushmail.com> wrote in message <hrd0jl$bsc$1@canopus.cc.umanitoba.ca>...
> ......
> Xc = (1/2) * a1 + (1/2) * a2,
> .....

  There is something decidedly fishy about those answers, Walter. In each of the answers 2, 3, and 4, they gave Xc as the average between a1 and a2. However, it would be easy to place the points (a1,b1) and (a2,b2) on an ellipse so that their x-coordinates are well off to one side of the ellipse's center and not equal to its average. What gives?

Roger Stafford

Subject: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)

From: Roger Stafford

Date: 29 Apr, 2010 23:07:04

Message: 8 of 18

  Gaurav, I should mention that the solution you obtain for the five coefficients may not actually represent an ellipse but rather a hyperbola or a parabola. If that is the case and there is no other column of S which has a zero singular value, then presumably there exists no ellipse which would run through those four given points. What you actually have here is the solution to a general four-point conic which is aligned with the coordinate axes in which the conic can be any one of these three types (or in some cases a degenerate conic.)

Roger Stafford

Subject: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)

From: GAURAV

Date: 29 Apr, 2010 23:26:04

Message: 9 of 18

Hi Roger and Walter,

In my case the points are at the extremes. So, two points on either side are at the border of the volume which is going to be formed.

               3
1 2 ... four points could be like this or
               4



    3
1 2... I doubt if some ellipse could be fit to these set of points..



     4

Anyways, I will try with the solved solution set provided by Walter and see what happens.
Any one of the last three solution set would work ?

And Roger, I need to form a surface enclosing it.. Be it ellipse or some closed surface. But on each slice along the height I have just four points to form a surface.
Any other ideas?

I wish I could give the small code to you and you could see the output and that would have made more sense. Can I have email id of both of you so that I can email you, if you guys dont mind?

Thanks a lot for your help.
Appreciate it.
Regards,
Gaurav.

Subject: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)

From: Bruno Luong

Date: 29 Apr, 2010 23:28:04

Message: 10 of 18

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hrd0f4$4o1$1@fred.mathworks.com>...
> "GAURAV " <gsharda@engineering.uiowa.edu> wrote in message <hrct9s$bpa$1@fred.mathworks.com>...
> > Hi Roger,
> > Thanks for your reply..
> >
> > Actually its an assumption I am making.
> > I have two curves(orthogonal to each other) aligned in 3D ( not necessarily symmetrical ).
> > The height (along Z axis) is the same. So, when I cut a slice for some Z, I have four points (two from each curve).
> > Now I need to construct a volume from these two views. So from each slice I am trying to fit in an ellipse. I hope this makes sense.
> >
> > So, now I have four points which pass through ellipse and ellipse should not go out of this bound .
> > So, I am looking to solve the ellipse equation and try to find out of the a,b,xc and yc required for the ellipse with angle=0.
> >
> > Can you suggest something?
> >
> > Regards,
> > Gaurav.
>
> For the four-point problem you can avoid using one of the non-linear methods by writing your equation in the following linear form. Then use 'svd' to determine the five unknown coefficients:
>
> A*x^2 + C*y^2 + D*x + E*y + F = 0
>
> Do this:
>
> x = [a1;a2;a3;a4];
> y = [b1;b2;b3;b4];
> M = [x.^2,y.^2,x,y,ones(4,1)];
> [U,S,V] = svd(M);
> A = V(1,5);
> C = V(2,5);
> D = V(3,5);
> E = V(4,5);
> F = V(5,5);
>

A note on Roger's solution: I would prefer normalize the eqt by fixing F = 1, for example. This will replace the calculation of svd by solving a 4x4 linear system.

Bruno

Subject: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)

From: Walter Roberson

Date: 29 Apr, 2010 23:42:38

Message: 11 of 18

Roger Stafford wrote:
> Walter Roberson <roberson@hushmail.com> wrote in message
> <hrd0jl$bsc$1@canopus.cc.umanitoba.ca>...
>> ......
>> Xc = (1/2) * a1 + (1/2) * a2,
>> .....
>
> There is something decidedly fishy about those answers, Walter. In
> each of the answers 2, 3, and 4, they gave Xc as the average between a1
> and a2. However, it would be easy to place the points (a1,b1) and
> (a2,b2) on an ellipse so that their x-coordinates are well off to one
> side of the ellipse's center and not equal to its average. What gives?


Good question.

When I tried running the original forms through my symbolic package, it
was taking an indefinite time. I normalized the forms and although the
solution wasn't exactly the fastest, it didn't take unduely long. But it
could be that I miswrote the normalization. If I have time, I'll try
again later... but I have to finish my taxes by tomorrow and tomorrow is
also the day my license for the symbolic package expires.

Subject: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)

From: GAURAV

Date: 29 Apr, 2010 23:44:05

Message: 12 of 18

Hi Walter and Roger,

Hey guys! The solution of four equations you gave works MUCH MUCH BETTER than the earlier ellipse fit thing I was using.

Thanks a lot!!! REALLY REALLY APPRECIATE.. Still it is not the best fit obtained but I gotta a hell lotta improvement.

Also, if something else you guys think would work for me to form a surface, let me know PLEASE!!!

Thank You,

Regards,
Gaurav.

Subject: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)

From: GAURAV

Date: 30 Apr, 2010 00:28:04

Message: 13 of 18

Hey Walter,

Can you tell me how you solved these four equations?

Thanks,
Gaurav

Subject: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)

From: Roger Stafford

Date: 30 Apr, 2010 00:28:04

Message: 14 of 18

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hrd4m4$2h1$1@fred.mathworks.com>...
> A note on Roger's solution: I would prefer normalize the eqt by fixing F = 1, for example. This will replace the calculation of svd by solving a 4x4 linear system.
>
> Bruno

  Hi Bruno. Normalizing by setting F to 0 doesn't always work. You can very easily be dealing with an ellipse that runs right through the origin, and the right value for F is zero in such a case. Even if it runs near the origin, forcing F to be 1 will encounter awkward accuracy problems.

  Using 'svd' has the virtue that the sum of the squares of the five coefficients is constrained to be unity which will easily allow F to be zero or some appropriately small value relative to the others. To me it seems the more "natural" method.

Roger Stafford

Subject: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)

From: Roger Stafford

Date: 30 Apr, 2010 00:59:04

Message: 15 of 18

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hrd86k$fje$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hrd4m4$2h1$1@fred.mathworks.com>...
> > A note on Roger's solution: I would prefer normalize the eqt by fixing F = 1, for example. This will replace the calculation of svd by solving a 4x4 linear system.
> >
> > Bruno
>
> Hi Bruno. Normalizing by setting F to 0 doesn't always work. You can very easily be dealing with an ellipse that runs right through the origin, and the right value for F is zero in such a case. Even if it runs near the origin, forcing F to be 1 will encounter awkward accuracy problems.
>
> Using 'svd' has the virtue that the sum of the squares of the five coefficients is constrained to be unity which will easily allow F to be zero or some appropriately small value relative to the others. To me it seems the more "natural" method.
>
> Roger Stafford

  Bruno, to make it seem a little less complicated, we could just as well have applied the 'null' function to the 4 x 5 M matrix I defined. We are simply looking for a non-zero vector which is orthogonal to its four rows in R^5, which would therefore satisfy the equation. There is certain to be one. Any one or more of the coefficients might naturally be zero without disrupting the calculations.

Roger Stafford

Subject: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)

From: Walter Roberson

Date: 30 Apr, 2010 01:53:36

Message: 16 of 18

GAURAV wrote:
> Hey Walter,
>
> Can you tell me how you solved these four equations?

Using a symbolic package not directly supported by Matlab:


F := (X,Y) -> (X-Xc)^2/(a^2) + (Y-Yc)^2/(b^2) - 1;

#that's the basic equation, but it is more efficient to rationalize it
#and then get rid of the denominator

F1 := (X,Y) -> numer(normal(F(X,Y));

solve([F1(a1,b1),F1(a2,b2),F1(a3,b3),F1(a4,b4)],[a,b,Xc,Yc]);

Based upon the length of time it is taking, I suspect that when I did
this earlier, I must have mis-written the equations. I'll post later
when it comes up with an answer.

Subject: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)

From: Bruno Luong

Date: 30 Apr, 2010 06:35:14

Message: 17 of 18

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hrd86k$fje$1@fred.mathworks.com>...

>
> Even if it runs near the origin, forcing F to be 1 will encounter awkward accuracy problems.

Roger, I don't think Roger it should not be any accuracy problem by fixing F=1. It is like divide the original implicit ellipsoid equation by F, and the linear equation to be solved is:

 M = [x.^2,y.^2,x,y];
X = - M(:,1:4) \ ones(4,1)

A = X(1)
C = X(2)
D = X(3)
E = X(4)
F = 1

In contrary, it should be more accurate (if not equal), because the new matrix above is submatrix of the 5-column matrix, thus the ratio of singular values becomes smaller because it is the projection.

I have my full code to fit ellipsoid in ndimensional that I have posted long ago in the newsgroup ( http://www.mathworks.com/matlabcentral/newsreader/view_thread/168603 ) where I just revised for axis-alignement:

function [H X0 W] = solveellipse(X, axisaligned)
% function [H X0 W] = solveellipse(X)
%
% function [H X0 W] = solveellipse(X, axisaligned)
%
% INPUT:
% X : (m x n) : m data points in R^n, they are supposed to belong
% to an ellipsoide (or more generally a second-order
% implicite hyper-surface)
% if axisaligned is TRUE, elleipsoid is forced to be aligned with axis
%
% OUTPUT:
% H : Matrix of the Bilinear form associated with the ellipsoide
% Elippsoide = { X in R^n : (X-X0)' * H * (X-X0) = 1 }
% X0 : (n x 1), center of the ellipsoide
% W : (n x n) square matrix where each column is the axis of ellipsoide
%

[ndata ndim]=size(X);

if ndim<2
    error('solveellipse: dimension number must be greater than 1');
end

% This vector will be used in few places for reshapping purpose
uno = ones(1,ndim);

%
% Generate all combinations of polynomial-order <=2 for ndim variables
%
order=cell(1,ndim);
order(:)={(0:2)};
ORDER=cell(1,ndim); % 1 x ndim cell of vectors (0:2)
[ORDER{:}]=ndgrid(order{:}); % Set {(0:2)}^ndim
ORDER=reshape(cat(ndim+1,ORDER{:}),[],ndim);
ORDER=ORDER(sum(ORDER,2)<=2,:); % second order only

%
% Remove constant term
%
ORDER=ORDER(2:end,:); % similar to ORDER(~any(ORDER,2),:)=[];

% Remove cross-term
if nargin>=2 && axisaligned
    ORDER(sum(ORDER==1,2)==2,:) = [];
end

npol=size(ORDER,1); % number of basis 2nd order-polynomials


if ndata<npol
    error('solveellipse: Not enough data');
end

    function loc=findloc(order) % nested function
        %
        % Look for the location of the 'order' in the set of {basis
        % 2nd-order polynomials}
        %
        [~, loc]=ismember(reshape(order,1,[]),ORDER,'rows');
    end

%
% Both of these are 3d array of (ndata x ndim x npol)
%
BIGX=repmat(X,[1 1 npol]);
BIGORDER=repmat(reshape(ORDER',[1 ndim npol]),[ndata 1 1]);
M=squeeze(prod(BIGX.^BIGORDER,2)); % product of power in every dimensions
clear BIGX BIGORDER

%
% Solve for polynomial coefficiens that lead to 1 on input points
%
rhs=ones(ndata,1);
P=M\rhs;
clear M rhs;

%
% Extract the Hessian from the solution
%
I=repmat(eye(ndim),ndim,1);
I=reshape(I,ndim*[1 1 1]);
J=permute(I,[2 1 3]); % swap the first two dimensions
K=squeeze(mat2cell(I+J,uno,uno,ndim));
P2_loc=cellfun(@findloc,K);
P2 = zeros(size(P2_loc));
P2(P2_loc>0) = P(P2_loc>0); % second order term
A=(tril(P2)+triu(P2))/2; % devide by 2, accept diagonal terms

%
% Extract the gradient
%
I=mat2cell(eye(ndim),uno,ndim);
P1_loc=cellfun(@findloc,I);
P1=-0.5*P(P1_loc); % -0.5 * first order term

X0=A\P1;
lambda=1/(1+X0'*P1);
H=lambda*A;

if nargout>=3 % Compute main-axis of ellipsoide
    [V D]=eig(H);
    d=diag(D);
    if any(d<0)
        warning('solveellipse:npdH', ...
                'solveellipse: non positive Hessian matrix');
    end
    W = V.*repmat(1./sqrt(d'),ndim,1);
end

end

% Bruno

Subject: SOLVE 4 NONLINEAR EQUATIONS (HELP PLEASE)

From: Roger Stafford

Date: 30 Apr, 2010 23:16:04

Message: 18 of 18

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hrdtn2$iit$1@fred.mathworks.com>...
> Roger, I don't think Roger it should not be any accuracy problem by fixing F=1. It is like divide the original implicit ellipsoid equation by F, and the linear equation to be solved is:
>
> M = [x.^2,y.^2,x,y];
> X = - M(:,1:4) \ ones(4,1)
>
> A = X(1)
> C = X(2)
> D = X(3)
> E = X(4)
> F = 1
>
> In contrary, it should be more accurate (if not equal), because the new matrix above is submatrix of the 5-column matrix, thus the ratio of singular values becomes smaller because it is the projection.
> ........
- - - - - - - - - -
  Hello Bruno. I ran a test of many randomly generated ellipses determined by four points so as to pass through the origin, and indeed what you stated seems to be true. In such situations, the

 "-[x.^2,y.^2,x,y]\ones(4,1)"

method appears to do as well as the

 "null([x.^2,y.^2,x,y,ones(4,1)])"

method on the average, even though I sometimes spaced the four points rather closely together.

  That is surprising to me because, with the theoretical ellipse running precisely through the origin, in principle the value of F should be exactly zero. The four-element method does indeed produce extremely large values for the remaining coefficients in an effort to minimize the effect of that F value of one and is apparently successful in doing so. Nevertheless I remain uneasy using such an algorithm, depending as it does on wayward round off errors for not actually going precisely through the origin. If it did, that method would encounter a truly singular matrix whose results in principle should have infinite values.

  As it is on my machine, with such ellipses that method always produces the disturbing notice: "Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = [some exceedingly small number]." It is easy to believe that the message might eventually be right in some special circumstances which I haven't managed to test.

  I suspect if and when the same kind of problem arises again I shall be sorely tempted to continuing using the 'null' function to find solutions for this homogenous type of equation. It seems more appropriate to my mathematical "soul".

Roger Stafford

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us