matlab怎么拟合出椭圆的函数啊 求解啊 急急急
数据拟合的最小二乘方法源于天文学中对行星或慧星这类天体的轨道计算。1795年,高斯在计算行星的椭圆轨道时提出并使用了这种方法,这一方法由勒让德于1805年首次公布。由开普列的研究成果,行星在其轨道平面上的运行轨迹是一个椭圆,而椭圆方程a1x2 + 2a2 xy + a3y2 + a4x + a5y + 1=0
需要由五个参数确定。原则上只要对行星的位置作5次观测就足以确定它的整个轨迹方程。但由于存在测量误差,由5次观测所确定的轨迹极不可靠,需要进行多次观测,用最小二乘法来消除误差,得到有关轨迹参数的更精确的值。最小二乘近似将几十次甚至上百次的观察所产生的高维空间问题降维到五参数的椭圆轨迹模型的五维空间处理。行星位置的10个观测点数据如下
x 1.02 0.95 0.87 0.77 0.67 0.56 0.44 0.30 0.16 0.01
y 0.39 0.32 0.27 0.22 0.18 0.15 0.13 0.12 0.13 0.15
利用最小二乘法确定行星的轨道参数,并确定椭圆方程。 这个应该不难,它形式上是关于x和y的二次函数,但是当把各观测值带入进去后,它是关于a1 a2 a3 a4 a5的一组线性方程,利用最小二乘原理即可解除其最小二乘解。 liwenhui 发表于 2012-8-22 10:25 static/image/common/back.gif
这个应该不难,它形式上是关于x和y的二次函数,但是当把各观测值带入进去后,它是关于a1 a2 a3 a4 a5的一组 ...
嗯,我也考虑用最小二乘法,百度找了好多程序,但是程序老是出错····能不能帮我看下 本帖最后由 liwenhui 于 2012-8-22 16:17 编辑
X.w.j.拽. 发表于 2012-8-22 10:39 static/image/common/back.gif
嗯,我也考虑用最小二乘法,百度找了好多程序,但是程序老是出错····能不能帮我看下
你需要首先明白最小二乘解的理论基础,然后才能动手写程序,它就是一个简单的矩阵运算:
function =stltp(x)
%- X,a Nx2 matrix,is the vector of satellite's position
%- PAR,a 2x2 matrix,is the parameters vector of satellite's position eqution
%- GM is the coeffient matrix of eqution
%- STLTP ueses least square method to estimate the the parameters vector of satellite's position eqution
=size(x);
if lx < 5,disp('input data is not enough'),return,end
if ly ~= 2,disp('input data dimension is not available'),return,end
xsq=x(:,1).^2;
xym=2*x(:,1).*x(:,2);
ysq=x(:,2).^2;
xp=x(:,1);
yp=x(:,2);
gm=;
b=-ones(lx,1);
gms=inv(gm'*gm);
par=gms*gm'*b;
调用这个函数: =stltp(mj')
得到:
par =
2.2538
0.0032
5.5222
-1.2898
-7.3774
gm =
1.0404 0.7956 0.1521 1.0200 0.3900
0.9025 0.6080 0.1024 0.9500 0.3200
0.7569 0.4698 0.0729 0.8700 0.2700
0.5929 0.3388 0.0484 0.7700 0.2200
0.4489 0.2412 0.0324 0.6700 0.1800
0.3136 0.1680 0.0225 0.5600 0.1500
0.1936 0.1144 0.0169 0.4400 0.1300
0.0900 0.0720 0.0144 0.3000 0.1200
0.0256 0.0416 0.0169 0.1600 0.1300
0.0001 0.0030 0.0225 0.0100 0.1500
其中向量par就对应了a1 a2 a3 a4 a5的值
我没有详细地验证结果对不对,不过思路应该没问题。
-------------------------------------------------------------------------------------------------------------------------
我刚才看了一下拟合后的效果,还不错,数据点都几乎分布在椭圆轨迹上:
{:3_55:}我回复的帖子哪去了? liwenhui 发表于 2012-8-22 19:07 static/image/common/back.gif
我回复的帖子哪去了?
贴子有在呀,很感谢你的回复。追问一下,如果X Y 的值增加到180个 那么程序改怎么改呢 X.w.j.拽. 发表于 2012-8-22 10:39 static/image/common/back.gif
嗯,我也考虑用最小二乘法,百度找了好多程序,但是程序老是出错····能不能帮我看下
不错,学习了!! 本帖最后由 liwenhui 于 2012-8-23 09:17 编辑
X.w.j.拽. 发表于 2012-8-22 22:15 static/image/common/back.gif
贴子有在呀,很感谢你的回复。追问一下,如果X Y 的值增加到180个 那么程序改怎么改呢
程序不用改,我提供的是一个m函数,你把m文件放在matlab搜索路径某个文件夹里,然后直接调用stltp(mj)就行,mj是N×2的一个矩阵,它的第一列表示x,第二列表示y。 liwenhui 发表于 2012-8-23 09:13 static/image/common/back.gif
程序不用改,我提供的是一个m函数,你把m文件放在matlab搜索路径某个文件夹里,然后直接调用stltp(mj)就 ...
我理解了你的意思,不过运行的时候我输入=stltp(mj')
mj=[2762 1792
2754 1786
2745 1781
2735 1775
2723 1769
2710 1764
2695 1758
2679 1753
2661 1747
2643 1742
2623 1737】;
出现这个错误??? Undefined function or variable 'mj'.
X.w.j.拽. 发表于 2012-8-23 09:57 static/image/common/back.gif
我理解了你的意思,不过运行的时候我输入=stltp(mj')
mj=[2762 1792
2754 1786
怎么可能用不了,你定义的mj最后的中括号用错了。
mj=[2762 1792
2754 1786
2745 1781
2735 1775
2723 1769
2710 1764
2695 1758
2679 1753
2661 1747
2643 1742
2623 1737]
=stltp(mj)
结果:
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =
1.610246e-19.
> In stltp at 16
a =
0.0000
-0.0000
0.0000
0.0002
-0.0015
b =
7628644 9899008 3211264 2762 1792
7584516 9837288 3189796 2754 1786
7535025 9777690 3171961 2745 1781
7480225 9709250 3150625 2735 1775
7414729 9633974 3129361 2723 1769
7344100 9560880 3111696 2710 1764
7263025 9475620 3090564 2695 1758
7177041 9392574 3073009 2679 1753
7080921 9297534 3052009 2661 1747
6985449 9208212 3034564 2643 1742
6880129 9112302 3017169 2623 1737
你的数据算出来,线性方程组的系数矩阵的转置对称逆矩阵是奇异矩阵,方程组的解接近奇异解,也就是系数矩阵之间存在高度共线性,最好换一组数据,不然结果无意义。
页:
[1]
2