数学建模社区-数学中国

标题: matlab怎么拟合出椭圆的函数啊求解啊急急急 [打印本页]

作者: X.w.j.拽. 时间: 2012-8-21 21:43
标题: 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
利用最小二乘法确定行星的轨道参数，并确定椭圆方程。

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

作者: X.w.j.拽. 时间: 2012-8-22 10:39

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

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

作者: liwenhui 时间: 2012-8-22 13:02
本帖最后由 liwenhui 于 2012-8-22 16:17 编辑

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 点体力)

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

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

作者: liwenhui 时间: 2012-8-22 19:07

我回复的帖子哪去了？

作者: X.w.j.拽. 时间: 2012-8-22 22:15

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

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

作者: wish_豪 时间: 2012-8-23 08:50

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

不错，学习了！！

作者: liwenhui 时间: 2012-8-23 09:13
本帖最后由 liwenhui 于 2012-8-23 09:17 编辑

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

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

作者: X.w.j.拽. 时间: 2012-8-23 09:57

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 时间: 2012-8-23 10:20

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

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

作者: X.w.j.拽. 时间: 2012-8-23 10:39

liwenhui 发表于 2012-8-23 10:20
怎么可能用不了，你定义的mj最后的中括号用错了。

mj=[2762 1792

听你这么说，我理解得七七八八了，不过，a1x2 + 2a2 xy + a3y2 + a4x + a5y + 1=0如果改为a1x2 + 2a2 xy + a3y2 + a4x + a5y + a6=0，那m文件了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];该怎么改呢，我改了然后运行不出，我们今天交作业，急得要死，真是麻烦你了。

作者: liwenhui 时间: 2012-8-23 11:41

X.w.j.拽. 发表于 2012-8-23 10:39
听你这么说，我理解得七七八八了，不过，a1x2 + 2a2 xy + a3y2 + a4x + a5y + 1=0如果改为a1x2 + 2a2 xy ...

这种情况是完全没必要考虑的，因为此时在方程两端同除以a6就化为和之前的模型完全一样了。也就是说，a1x2 + 2a2 xy + a3y2 + a4x + a5y + 1=0和a1x2 + 2a2 xy + a3y2 + a4x + a5y + a6=0是一个模型，估计出来的参数a6肯定为1。

作者: X.w.j.拽. 时间: 2012-8-23 12:22

liwenhui 发表于 2012-8-23 11:41
这种情况是完全没必要考虑的，因为此时在方程两端同除以a6就化为和之前的模型完全一样了。也就是说，a1 ...

嗯，解决了，谢谢你，真是高手

作者: 琪_He 时间: 2013-5-20 21:39

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

如果加上‘椭圆的焦点在原点’这一条件，该怎么处理呢

作者: liwenhui 时间: 2013-5-21 13:40

琪_He 发表于 2013-5-20 21:39
如果加上‘椭圆的焦点在原点’这一条件，该怎么处理呢

当焦点的坐标被确认后，代数上看实质是对最小二乘方法施加了一个线性约束条件，就可采用有约束的最小二乘法。

QQ截图20130521132420.jpg (14.32 KB, 下载次数: 201)

QQ截图20130521132420.jpg

作者: 琪_He 时间: 2013-5-24 20:28

liwenhui 发表于 2013-5-21 13:40
当焦点的坐标被确认后，代数上看实质是对最小二乘方法施加了一个线性约束条件，就可采用有约束的最小二乘 ...

椭圆 X坐标 5.764 6.286 6.759 7.168 7.480 Y坐标 0.648 1.202 1.832 2.526 3.360 并且焦点在原点，请问您怎么用最小二乘法条件拟合呢？

作者: Fey1992 时间: 2013-5-25 23:41

liwenhui 发表于 2012-8-22 13:02
你需要首先明白最小二乘解的理论基础，然后才能动手写程序，它就是一个简单的矩阵运算：

function = ...

大神手里有个交通流的元胞自动机自己水平低了点有点看不懂？可以指教下嘛

作者: jonejack11 时间: 2013-6-6 00:55
有学到了新知识

作者: 沐之城 时间: 2013-8-18 22:32
支持一下！

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