Thread Subject:

fitting curve at polar coordinate

I have some empirical data point at x-y coordinate from my research and a suggested equation in polar coordinate derived from a modeling. The equation after simplification is like (R=a/(1 + b cos(theta) ), where a and b are constants to be found. I want to evaluate a, b coefficients from my data.
My question is: Can I use MATLAB for this problem? If yes, how?

"M Srj" wrote in message <jigqp3$f6e$1@newscl01ah.mathworks.com>...
> I have some empirical data point at x-y coordinate from my research and a suggested equation in polar coordinate derived from a modeling. The equation after simplification is like (R=a/(1 + b cos(theta) ), where a and b are constants to be found. I want to evaluate a, b coefficients from my data.
> My question is: Can I use MATLAB for this problem? If yes, how?
Yes, this is a fairly simple problem to use Matlab to solve. To get a detailed description of what to do you should search the function "lsqcurvefit" in the help menu.
http://www.mathworks.com/help/toolbox/optim/ug/lsqcurvefit.html
You want to write a function where your input is theta and param where param = [a,b]
for example write this text in an m file and save it.
function R = fit(theta, param)
a = param(1);
b = param(2);
R=a/(1 + b cos(theta);
end

create vectors thetadata and Rdata
then in the command window type param = lsqcurvefit(@fit, [1 1], thetadata, Rdata)
this should return an leastsquares fit of a and b

Thanks Michael for your reply. I run this m-file but got error messages in MATLAB after my data in (theta, R) format:
 
Error using /
Matrix dimensions must agree.
Error in fit (line 4)
R=a/(1 + b* cos(theta*pi/180));
Error in lsqcurvefit (line 205)
            initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});
Caused by:
    Failure in initial user-supplied objective function evaluation. LSQCURVEFIT cannot continue.

Do know what wrong is?
Thanks
 
 
"Michael W." wrote in message <jki2uc$ocm$1@newscl01ah.mathworks.com>...
> "M Srj" wrote in message <jigqp3$f6e$1@newscl01ah.mathworks.com>...
> > I have some empirical data point at x-y coordinate from my research and a suggested equation in polar coordinate derived from a modeling. The equation after simplification is like (R=a/(1 + b cos(theta) ), where a and b are constants to be found. I want to evaluate a, b coefficients from my data.
> > My question is: Can I use MATLAB for this problem? If yes, how?
> Yes, this is a fairly simple problem to use Matlab to solve. To get a detailed description of what to do you should search the function "lsqcurvefit" in the help menu.
> http://www.mathworks.com/help/toolbox/optim/ug/lsqcurvefit.html
> You want to write a function where your input is theta and param where param = [a,b]
> for example write this text in an m file and save it.
> function R = fit(theta, param)
> a = param(1);
> b = param(2);
> R=a/(1 + b cos(theta);
> end
>
> create vectors thetadata and Rdata
> then in the command window type param = lsqcurvefit(@fit, [1 1], thetadata, Rdata)
> this should return an leastsquares fit of a and b

"M Srj" wrote in message <jn22c5$bab$1@newscl01ah.mathworks.com>...
> Thanks Michael for your reply. I run this m-file but got error messages in MATLAB after my data in (theta, R) format:
> .......
> Error using /
> Matrix dimensions must agree.
> Error in fit (line 4)
> R=a/(1 + b* cos(theta*pi/180));
> ......
> "Michael W." wrote in message <jki2uc$ocm$1@newscl01ah.mathworks.com>...
> > function R = fit(theta, param)
> > a = param(1);
> > b = param(2);
> > R=a/(1 + b cos(theta);
> > end
- - - - - - - - - -
  The function you pass to 'lsqcurvefit' should have the parameters as its first, not its second, argument: R = fit(param,theta). Also unless I am greatly mistaken, 'lsqcurvefit' passes 'fit' a vector of theta values, in which case you should have:

 R=a./(1 + b cos(theta);

with the "./" kind of division in the 'fit' definition.

  The curve you are trying to fit your data to is an ellipse, so as an alternative you might want to check out some of the numerous threads in this newsgroup that have dealt with fitting ellipses. You might get a far better fit that way.

Roger Stafford

"Roger Stafford" wrote in message <jn27iu$h0$1@newscl01ah.mathworks.com>...
> The curve you are trying to fit your data to is an ellipse, .....
- - - - - - - - -
  Correction. I should have said your curve model is either an ellipse or a hyperbola depending on whether b^2 is less than or greater than 1, respectively.

Roger Stafford

Thanks! It is working now and give me parameters a and b.
However, the curve fitting is very poor! Actually, my experimental data are like parabola. Do you know better fitting way?

"M Srj" wrote in message <jn4dlf$i0a$1@newscl01ah.mathworks.com>...
> Thanks! It is working now and give me parameters a and b.
> However, the curve fitting is very poor! Actually, my experimental data are like parabola. Do you know better fitting way?
- - - - - - - - - - -
  It's not surprising that your fit is not good, because you have allowed yourself only two degrees of freedom in the family of functions to fit with R = a/(1+b*cos(theta)). If you allow the full five parameters of general conics, you could probably get a much better fit. Here is an approach you might consider.

  Any conic section (including parabolas) can be characterized by the following equation:

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

where A, B, C, D, E, and F are appropriate constants. Let x and y be column vectors of the pairs of cartesian coordinates of the data you wish to fit. Then do this:

 M = [x.^2,x.*y,y.^2,x,y,ones(size(x))];
 [U,S,V] = svd(M,0);

The sixth (last) column of V will contain the least squares values of [A;B;C;D;E;F] for the above expression subject to the restriction that the sum of their squares is one. You can test how good the fit is by computing

 A*x.^2 + B*x.*y + C*y.^2 + D*x + E*y + F

to see how near zero it is for your data.

  If your data really is an approximate parabola, you should find that B^2 will be approximately equal to 4*A*C.

  Of course this still leaves you with the problem of how to generate a curve from the above equation, but I'll leave you with that problem for now.

Roger Stafford

I used svd function and found U, S, and V. I looked at MATLAB help for "singular value decomposition" but found no clue for finding A, B, ,,,F!

"M Srj" wrote in message <jn4slf$h8f$1@newscl01ah.mathworks.com>...
> I used svd function and found U, S, and V. I looked at MATLAB help for "singular value decomposition" but found no clue for finding A, B, ,,,F!
- - - - - - - - -
  They are in the sixth column of V:

 A = V(1,6)
 B = V(2,6)
 ...
 F = V(6,6)

  This is the eigenvector possessing the smallest singular value, which in turn means the expression A*x.^2+B*x.*y+ ... is closest to zero.

Roger Stafford

Okey!
Now, I have polynomial A*x.^2+B*x.*y+ .... which is almost zero at my X, Y data.
Two questions
1- What is the fitted curve to my data?
2- How to plot my data and the fitted curve in one figure?
Thanks

"M Srj" wrote in message <jn52gu$9vo$1@newscl01ah.mathworks.com>...
> Okey!
> Now, I have polynomial A*x.^2+B*x.*y+ .... which is almost zero at my X, Y data.
> Two questions
> 1- What is the fitted curve to my data?
> 2- How to plot my data and the fitted curve in one figure?
> Thanks
- - - - - - - - -
  (I was afraid you would ask me that!) The equation

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

is itself an implicit form of the fitted curve which you are seeking. The most direct way to generate the curve for plotting purposes is to solve either x or y in terms of the other. This will in general lead to solving a quadratic equation which will have either two solutions, no solutions, or at special points one solution.

  As an example, let's suppose you decide to solve for y in terms of x in the equation

 5*x^2-6*x*y+2*y^2-x+5*y+3=0

Rewrite it as:

 2*y^2 + (-6*x+5)*y + (5*x^2-x+3) = 0

By the well-known quadratic equation formula the two solutions would then be:

 y1 = (6*x-5+sqrt((6*x-5)^2-4*2*(5*x^2-x+3)))/4 and
 y2 = (6*x-5-sqrt((6*x-5)^2-4*2*(5*x^2-x+3)))/4

The quantity inside the square root can be expressed as:

 (6*x-5)^2-4*2*(5*x^2-x+3)) = -4*x^2-62^x=1 = -4*(x+13/2)^2+170

which means that

 (x+13/2)^2 <= 170/4
 -1/2*sqrt(170)-13/2 <= x <= 1/2*sqrt(170)-13/2

In other words what we have here is an ellipse and x must range between
these two numbers. We need to plot both y1 and y2 between these two limits to get the entire ellipse. (Use X for this x to distinguish it from the original data x.)

 X = linspace(-1/2*sqrt(170)-13/2+1e-14,1/2*sqrt(170)-13/2-1e-14,500);
 y1 = (6*X-5+sqrt((6*X-5).^2-4*2*(5*X.^2-X+3)))/4;
 y2 = (6*X-5-sqrt((6*X-5).^2-4*2*(5*X.^2-X+3)))/4;
 plot(X,y1,'y-',X,y2,'y-',x,y,'ro') % Plot the ellipse and the data points
 axis equal

  There are other ways of generating conic section curves which don't require two distinct parts as above but rather use a parameter. I would prefer not going into that now, however.

Roger Stafford

Thanks for your helps.
I found a general form for downward parabola in polar coordinate as R=a.b/(1+b.sin(theta))
and calculated parameters a and b by using "lsqcurvefit" command and the m.file .
Now, the fitting is very good in polar form.

Correction: equation R=a.b/(1+b.sin(theta)) is a general form for conic section by vertical axis of symmetry not only opening downward parabola.

"M Srj" wrote in message <jn6ssl$c8q$1@newscl01ah.mathworks.com>...
> I found a general form for downward parabola in polar coordinate as R=a.b/(1+b.sin(theta))
- - - - - - - -
  That is by no means a general form for parabolas. To be a parabola, you must have b^2 equal to 1 which leaves only one degree of freedom. General parabolas have four. "Downward" parabolas have three.

Roger Stafford

I already said: "Correction: equation R=a.b/(1+b.sin(theta)) is a general form for conic section by vertical axis of symmetry not only opening downward parabola." which was posted before your reply. You can check this link for above equation:
http://jwilson.coe.uga.edu/EMAT6680Fa07/Cody/Assignment11/Assignment11.html

I do know parameter b=1 for parabola but the above equation fits well to my experimental data with two parameters a and b not equal one. So, my data may be a part of an ellipse (I just realized it!).
Anyway, do you have the equation with more DOF in Polar coordinate or Cartesian coordinate? Because my main goal is to develop a correlation out of my data in terms of R and theta.
Thanks

"M Srj" wrote in message <jn711v$2ho$1@newscl01ah.mathworks.com>...
> I already said: "Correction: equation R=a.b/(1+b.sin(theta)) is a general form for conic section by vertical axis of symmetry not only opening downward parabola." which was posted before your reply. You can check this link for above equation:
> http://jwilson.coe.uga.edu/EMAT6680Fa07/Cody/Assignment11/Assignment11.html
>
> I do know parameter b=1 for parabola but the above equation fits well to my experimental data with two parameters a and b not equal one. So, my data may be a part of an ellipse (I just realized it!).
- - - - - - - - - - -
  A curve with the polar equation r = a*b/(1+b*sin(theta)) is a conic section with the origin of the polar coordinates located at one of the focal points of the conic, with eccentricity 'b', and with the directrix parallel to the x-axis and a distance 'a' above it. It is not a general conic since a conic might have its directrix in any orientation in a full circle and with a focal point located elsewhere than at the origin.

  Thus a set of points you are trying to fit might have its best fitting conic with focal point away from the origin and directrix oriented in some other direction than parallel to the x-axis. If you actually found a good fit for your data with the above equation, that means you were fortunate.

Roger Stafford

Thanks for all your information. I learnt a lot from your answers.
Actually, I studied the conditions to this equation and then chose it. My data present downward open curve smooth contour with vertical symmetry.
Anyway, you said of downward parabola equation with 3 DOF. But my equation for parabola has only 1 DOF. Is your equation with 3 DOF in Polar or Cartesian coordinate? and what is it?
 M

"M Srj" wrote in message <jn7fua$2fm$1@newscl01ah.mathworks.com>...
> Thanks for all your information. I learnt a lot from your answers.
> Actually, I studied the conditions to this equation and then chose it. My data present downward open curve smooth contour with vertical symmetry.
> Anyway, you said of downward parabola equation with 3 DOF. But my equation for parabola has only 1 DOF. Is your equation with 3 DOF in Polar or Cartesian coordinate? and what is it?
> M
- - - - - - - - - - -
  I don't have the time to give detailed equations or matlab code as an answer to your question so I will describe in words the five-parameter case that would lead to a polar form of a fitted curve. It would produce a function that you could presumably pass to 'lsqcurvefit' for obtaining the best set of parameters, though admittedly I have not done this myself. (I have never done a search on the file exchange for this - maybe someone has already carried all this out.)

  Starting with data points in cartesian coordinates, a pair of parameters, d and e, would be the x and y locations of a presumed focal point of a conic. These would be subtracted from the data coordinates to accomplish a translation of origin. The next step would be conversion of these translated cartesian coordinates to polar coordinates, r and theta, for each data point. Finally these would be placed in the expression

 r*(1+b*sin(theta-c))-a*b

or if you prefer

 r*(1+b*cos(theta-c))-a*b

containing three further parameters, a, b, and c (note the addition of parameter c here.) These five parameters, a, b, c, d, and e, would be adjusted by 'lsqcurvefit' in order to obtain a least squares average of the above expression over all the points. The five parameters are a: distance to directrix, b: eccentricity, c: angle of orientation of the directrix, d & e: cartesian coordinates of the conic focal point.

  When you finish, if successful, you would have the five best parametric values and these could then be used to generate a plot of the conic by using theta as a parameter varied through a 2*pi range, solving for r in the conic equation, and then transforming back to the original cartesian coordinates for plotting purposes.

  I apologize for the vagueness in the above description. I can visualize it all but unfortunately don't have time to write it all out for you.

Roger Stafford