本帖最后由 xuhai777 于 2018-10-13 21:42 编辑
已知N个样本点,用多项式拟合过这N个点的曲线.
理论的话就不多说了.
n次多项式:Y=ax^n+bx^(n-1)...fx^1+g
已知n次多项式的(n+1)个样本,用消元法求解多项式的系数
思路:
1)多项式按X的幂由低到高排序格式
a+bx+cx^2...+gx^n=Y
2)转换格式:n次则有(n+1)系数,加上Y值,所以列数是(n+2)
a,b,c...g,y
3)把(n+1)个样本组数分别代入,并计算X的幂,书写矩阵格式,行数是(n+1)
4)从左到右顺序设置主系数,开始消元并归1化,直到第n系数先被计算出来
消元:已归1化后的行,异行同列相减,可使主系数为0
归1化:行内全部系数除以主系数
5)把n系数回代入(n-1)行,计算(n-1)系数
再把n和(n-1)系数回代入(n-2)行,计算(n-2)系数
依此法,直到计算出全部系数,系数在末列
限制条件:
样本唯一的x值对应唯一的y值,也就是多项式的曲线不能分叉
例举:已知3点{1,2 ,2,4 ,4,3},求过此3点的多项式曲线
解法:设此曲线的多项式为 a+bx+cx^2=Y
(1)分别代入3点数据,列出矩阵式
1 1 1 2 行1
1 2 4 4 行2
1 4 16 3 行3
(2)消元并归1化,计算出n系数
行2=行2-行1 ;行3=行3-行1
1 1 1 2 行1
=> 0 1 3 2 行2
0 3 15 1 行3
行3归1化
1 1 1 2 行1
=> 0 1 3 2 行2
0 1 5 1/3 行3
行3=行3-行2
1 1 1 2 行1
=> 0 1 3 2 行2
0 0 2 -5/3 行3
行3归1化
1 1 1 2 行1
=> 0 1 3 2 行2
0 0 1 -5/6 行3
(3)回代,直到计算出全部系数
行3代入行2
1 1 1 2 行1
=> 0 1 0 9/2 行2
0 0 1 -5/6 行3
行3行2代入行1
1 0 0 -5/3 行1
=> 0 1 0 9/2 行2
0 0 1 -5/6 行3
(4)所求多项式
-5/3 +9/2x -5/6x^2 =Y
C语言代码
//float s[]={1,2 ,2,4 ,4,3 ,5,2 ,7,-1 ,7.5,0.5}; //样本:x,y
float s[]={1,2 ,2,4 ,4,3 ,5,2 ,7,-1 ,6,0.5}; //样本:x,y
//float s[]={1,2 ,2,4 ,3,6}; //直线,计算结果中无用系数是0
#define n (sizeof(s)/sizeof(float)/2) //样本组数,行数
#define m (n+1) //矩阵行变量数
//---------------------------
float Pow(float i,int j) //幂函数,x^0=1;x^1=x
{ float k=1;
while(j--) { k*=i; };
return(k);
}
void poly(void)
{ int i,j,k,z,r;
float t[n*m],w; //n*(n+1),计算用的矩阵数组,计算结果在末列
if(n>1)
{ z=0; k=0; //初始化矩阵,顺序:x^0,x^1,x^2...y
for(j=0;j<n;j++)
{ w=s[k++]; //X值
for(i=0;i<n;i++) //第一列总是1,没优化
{ t[z++]=Pow(w,i); } //X的幂计算
t[z++]=s[k++]; //Y值
}
//-----------------------
for(k=0;k<n;k++) //总元数
{ z=k*m; //初行,参考行
for(j=k+1;j<n;j++) //行数
{ r=j*m; //始行,消元行
for(i=k;i<m;i++) //列数
{ t[r+i]-=t[z+i]; } //消元
r++;
if(t[r+k]!=1) //行元位置
{ for(i=k+1;i<n;i++) //列数
{ t[r+i]/=t[r+k]; } //行归1化
t[r+k]=1; //为了好看
}
}
}
//-----------------------
for(j=n-1;j>0;j--) //行
{ w=0;
for(i=j;i<n;i++) //列
{ w+=t[m*(j-1)+i]*t[m*(i+1)-1]; }
t[m*j-1]-=w; //计算行Y植
z=m*(j-1); for(i=j;i<n;i++) { t[z+i]=0; } //为了好看
}
}
while(1); //断点
}
//===========================
int main(void)
{ poly();
while(1);
} |
[复制链接]
发表于 2018-10-13 21:23:39
楼主
