matlab怎么拟合出椭圆的函数啊 求解啊 急急急

X.w.j.拽.

　　数据拟合的最小二乘方法源于天文学中对行星或慧星这类天体的轨道计算。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
利用最小二乘法确定行星的轨道参数，并确定椭圆方程。

liwenhui

这个应该不难，它形式上是关于x和y的二次函数，但是当把各观测值带入进去后，它是关于a1 a2 a3 a4 a5的一组线性方程，利用最小二乘原理即可解除其最小二乘解。

X.w.j.拽.

liwenhui 发表于 2012-8-22 10:25
这个应该不难，它形式上是关于x和y的二次函数，但是当把各观测值带入进去后，它是关于a1 a2 a3 a4 a5的一组 ...

嗯，我也考虑用最小二乘法，百度找了好多程序，但是程序老是出错····能不能帮我看下

liwenhui

X.w.j.拽. 发表于 2012-8-22 10:39
嗯，我也考虑用最小二乘法，百度找了好多程序，但是程序老是出错····能不能帮我看下

你需要首先明白最小二乘解的理论基础，然后才能动手写程序，它就是一个简单的矩阵运算：

function [par,gm]=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
    [lx,ly]=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=[xsq,xym,ysq,xp,yp];
    b=-ones(lx,1);
    gms=inv(gm'*gm);
    par=gms*gm'*b;

调用这个函数： [par,gm]=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的值

我没有详细地验证结果对不对，不过思路应该没问题。 stltp.m (665 Bytes, 下载次数: 2, 售价: 1 点体力)

2012-8-22 13:02 上传

点击文件名下载附件
售价: 1 点体力  [记录]
下载积分: 体力 -2 点

-------------------------------------------------------------------------------------------------------------------------

我刚才看了一下拟合后的效果，还不错，数据点都几乎分布在椭圆轨迹上：

liwenhui

我回复的帖子哪去了？

X.w.j.拽.

liwenhui 发表于 2012-8-22 19:07
我回复的帖子哪去了？

贴子有在呀，很感谢你的回复。追问一下，如果X Y 的值增加到180个 那么程序改怎么改呢

wish_豪

X.w.j.拽. 发表于 2012-8-22 10:39
嗯，我也考虑用最小二乘法，百度找了好多程序，但是程序老是出错····能不能帮我看下

不错，学习了！！

liwenhui

X.w.j.拽. 发表于 2012-8-22 22:15
贴子有在呀，很感谢你的回复。追问一下，如果X Y 的值增加到180个 那么程序改怎么改呢

程序不用改，我提供的是一个m函数，你把m文件放在matlab搜索路径某个文件夹里，然后直接调用stltp(mj)就行，mj是N×2的一个矩阵，它的第一列表示x，第二列表示y。

X.w.j.拽.

liwenhui 发表于 2012-8-23 09:13
程序不用改，我提供的是一个m函数，你把m文件放在matlab搜索路径某个文件夹里，然后直接调用stltp(mj)就 ...

我理解了你的意思，不过运行的时候我输入[par,gm]=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'.

liwenhui

X.w.j.拽. 发表于 2012-8-23 09:57
我理解了你的意思，不过运行的时候我输入=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]

[a,b]=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

你的数据算出来，线性方程组的系数矩阵的转置对称逆矩阵是奇异矩阵，方程组的解接近奇异解，也就是系数矩阵之间存在高度共线性，最好换一组数据，不然结果无意义。