有点的坐标求摆线的参数方程

laihai718

已知如下数据：
t=0  0.1   0.2   0.3   0.4   0.5   0.6   0.7    0.8   0.9   1
x =13.7836   13.9968   14.4546   15.1899   16.2206   17.5497   19.1641   21.0355   23.1201   25.3607   27.6875
y =12.9607   12.9333   13.0332   13.1918   13.3350   13.3855   13.2660   12.9017   12.2234   11.1702    9.6919
x=a*sin(b+c+t)-e*cos(f)*sin(b+c+g*t)
y=a*cos(b+c+t)-e*cos(f)*cos(b+c+g*t)
我用的是弧度制。
求x、y中的系数a,b,c,e,f,g。

请教方法，不胜感激

swinder

此貌似非线形拟和， lsqcurvefit还有nlinfit
个人比较喜欢后者。

shamohu

用1stOpt试试：

SharedModel;
Variable t,x,y;
Function
x=a*sin(b+c+t)-e*cos(f)*sin(b+c+g*t);
y=a*cos(b+c+t)-e*cos(f)*cos(b+c+g*t);
Data;
0        13.7836        12.9607
0.1        13.9968        12.9333
0.2        14.4546        13.0332
0.3        15.1899        13.1918
0.4        16.2206        13.3350
0.5        17.5497        13.3855
0.6        19.1641        13.2660
0.7        21.0355        12.9017
0.8        23.1201        12.2234
0.9        25.3607        11.1702
1        27.6875        9.6919

结果：
均方差(RMSE): 3.09420406252851E-5
残差平方和(RSS): 2.10630173172495E-8
相关系数(R): 0.999999999998214
相关系数之平方(R^2): 0.999999999996428

参数                  最佳估算
--------------------        -------------
a                  45.0003230755287
b                  -37.7876692351975
c                  38.6038273703678
e                  26.0835760471171
f                  -0.0158553828616697
g                  1.66665847505903

因为你公式的原因，上述结果除了a和g外，都不是唯一的。你的公式实际上可简化成：
x=a*sin(b+t)-c*sin(b+g*t);
y=a*cos(b+t)-c*cos(b+g*t);

laihai718

真实值是a =45，b =0.8378，c =-0.0216，e =27，f =0.2618，g = 1.6667

To:2楼swinder
我用您所提供的两种方法拟合，结果误差都很大，不知道是什么原因？

用 lsqcurvefit拟合：

function F = baixian(x,tdata,xdata)
F=x(1)*sin(x(2)+x(3)+tdata)-x(4)*cos(x(5))*sin(x(2)+x(3)+x(6)*tdata)-xdata;

tdata = [0;0.1;0.2;0.3;0.4;0.5;0.6;0.7;0.8;0.9;1];
xdata = [13.7836;13.9968;14.4546;15.1899;16.2206;17.5497;19.1641;21.0355;23.1201;25.3607;27.6875];
x = lsqcurvefit(@myfun, [40 0.8 0 25 0.05 1], tdata, xdata)

16.1905   11.4896         0   25.0000    0.0500    1.0000

用nlinfit拟合：
clear;
t=[0  0.1   0.2   0.3   0.4   0.5   0.6   0.7    0.8   0.9   1];
y=[13.7836   13.9968   14.4546   15.1899   16.2206   17.5497   19.1641   21.0355   23.1201   25.3607   27.6875];
myfunc=inline('x(1)*sin(x(2)+x(3)+t)-x(4)*cos(x(5))*sin(x(2)+x(3)+x(6)*t)','x','t');
x=nlinfit(t,y,myfunc,[0.5 0.5 0.5 0.5 0.5 0.5]);
a=x(1),b=x(2),c=x(3),e=x(4),f=x(5),g=x(6)

a =-24.2231，b =-1.3570，c =-1.3570，e =-3.5975，f =2.2026，g = -180.7152

laihai718

To：shamohu
您说的优化方法好像是个专门的优化软件，您说解出来有好几组答案，不知道怎么判别哪种是对哪个是错？

swinder

可以使用option，我原来要用上面的函数得到用最小二乘法拟和的解六位精度的解时，需要在option里设置很高的精度才达得到。你可以试哈

laihai718

按照shamohu版主的方法，我用LM+UGO算法，也是没运行一次，就出一个结果，a和g的结果几乎没变化，其他的都在变，下面这个结果是比较切近实际的结果：
a=44.9984296
b=0.6321948
c=0.18396917
e=26.0827369
f=0.014685
g=1.666671382

由于软件没有注册，不能复制

laihai718

另外，每次运行结果我还发现一个问题，其中a  e  g是实在的数，变化幅度不大，b  c  f是弧度，每次运行变化幅度都很大，会不会跟这个有关系？

laihai718

可以使用option，我原来要用上面的函数得到用最小二乘法拟和的解六位精度的解时，需要在option里设置很高的精度才达得到。你可以试哈
swinder 发表于 2009-10-13 14:13

请问option是个小软件？

swinder

不是，是上述函数的精度啊，跳出条件，步长啊等参数设置，help里有。

shamohu

另外，每次运行结果我还发现一个问题，其中a  e  g是实在的数，变化幅度不大，b  c  f是弧度，每次运行变化幅度都很大，会不会跟这个有关系？
laihai718 发表于 2009-10-13 15:00

三角函数是周期性的，不加约束，会有很多解。

laihai718

三角函数是周期性的，不加约束，会有很多解。
shamohu 发表于 2009-10-13 16:39

我加个约束试试

laihai718

请教下版主，如果我这个题没有t的值，该怎么解？

shamohu

如果t未知，可用如下代码：

SharedModel;
ParVariable t;
Variable x,y;
Function
x=a*sin(b+c+t)-e*cos(f)*sin(b+c+g*t);
y=a*cos(b+c+t)-e*cos(f)*cos(b+c+g*t);
Data;
13.7836        12.9607
13.9968        12.9333
14.4546        13.0332
15.1899        13.1918
16.2206        13.3350
17.5497        13.3855
19.1641        13.2660
21.0355        12.9017
23.1201        12.2234
25.3607        11.1702
27.6875        9.6919

hujianbiao403

可以用拟合工具箱进行拟合