数学建模社区-数学中国

标题: 在线MATLAB高手前来一看！！ [打印本页]

作者: 神帝在线 时间: 2009-8-6 13:38
标题: 在线MATLAB高手前来一看！！
x y
-32.303 18.082
-33.623 23.257
-37.209 26.523
-42.182 27.989
-47.117 26.457
-50.902 22.743
-52.23 18.26
-50.731 12.772
-47.117 26.457
-42.529 7.963
-36.87 9.408
-33.576 12.738
-32.496 19.968
-34.561 24.304
-38.462 27.116
-44.273 27.749
-48.935 25.409
-51.42 21.642
-52.09 16.517
-49.686 11.152
-45.783 8.528
-40.278 8.116
-35.54 10.533
-32.88 14.151
-32.622 21.015
-35.93 25.629
-40.631 27.703
-45.736 27.168
-49.97 24.216
-52.157 19.489
-51.594 14.296
-48.322 10.049
-44.249 8.184
-39.064 8.313
-34.761 11.227
-32.407 16.009

36个点的坐标，我们想用最小二乘法拟合成圆，求出圆心，并确定半径。编程如下：
function
[xc,yc,R,a] = circfit(x,y)
%
[xc yx R] = circfit(x,y)
%
圆心为 (yc,xc) 半径为 R
%
x^2+y^2+a*x+b*y+c=0
%
用最小二乘法，xc=-0.5a(1),yc=-0.5a(2)

x=x(; y=y(;

abc=[x y ones(size(x))]\[-(x.^2+y.^2)];

a=abc(1);

b=abc(2);

c=abc(3);

xc = -.5*a;

yc = -.5*b;

R
=
sqrt((a^2+b^2)/4-c)
th = linspace(0,2*pi,200)';
plot(x,y,'o'), title(' measured points')

pause(1)

% 由以上数据拟合圆

[xc,yc,Re,a] = circfit(x,y);

xe = Re*cos(th)+xc; ye = Re*sin(th)+yc;

plot(x,y,'o',[xe;xe(1)],[ye;ye(1)],'-.'),

title(' measured fitted and true circles')

legend('measured','fitted')

text(xc-Re*0.9,yc,sprintf('center (%g , %g );
R=%g',xc,yc,Re))

xlabel x, ylabel y

axis equa
我电脑的MATLAB打不开了，希望MATLAB高手运行下，出个结果。

欢迎光临数学建模社区-数学中国 (http://www.madio.net/)