多项式拟合

1. 已知n个坐标点，求出逼近曲线，就是说要用一个多项式曲线拟合这些离散点咯，有很多使用样条曲线进行拟合的
2. 得到曲线的多项式表达式，项数不限，这个问题就模糊了，项数越多，越逼近，如果是n个坐标点，用n项也就是n-1次多项式就能插值每个点了，系数通过求解一个线性方程组就可以了；不过，次数太高，不知意义何在

解决方案 »

免费领取超大流量手机卡，每月29元包185G流量+100分钟通话, 中国电信官方发货

最小二乘法曲线拟合：
typedef CArray<double,double>CDoubleArray;
BOOL CalculateCurveParameter(CDoubleArray *X,CDoubleArray *Y,long M,long N,CDoubleArray *A)
{
//X,Y --  X,Y两轴的坐标
//M   --  结果变量组数
//N   --  采样数目
//A   --  结果参数，从低次到高次

register long i,j,k;
double Z,D1,D2,C,P,G,Q;
CDoubleArray B,T,S;
B.SetSize(N);
T.SetSize(N);
S.SetSize(N);
if(M>N)M=N;
for(i=0;i<M;i++)
(*A)[i]=0;
Z=0;
B[0]=1;
D1=N;
P=0;
C=0;
for(i=0;i<N;i++)
{
P=P+(*X)[i]-Z;
C=C+(*Y)[i];
}
C=C/D1;
P=P/D1;
(*A)[0]=C*B[0];
if(M>1)
{
T[1]=1;
T[0]=-P;
D2=0;
C=0;
G=0;
for(i=0;i<N;i++)
{
Q=(*X)[i]-Z-P;
D2=D2+Q*Q;
C=(*Y)[i]*Q+C;
G=((*X)[i]-Z)*Q*Q+G;
}
C=C/D2;
P=G/D2;
Q=D2/D1;
D1=D2;
(*A)[1]=C*T[1];
(*A)[0]=C*T[0]+(*A)[0];
}
for(j=2;j<M;j++)
{
S[j]=T[j-1];
S[j-1]=-P*T[j-1]+T[j-2];
if(j>=3)
{
for(k=j-2;k>=1;k--)
S[k]=-P*T[k]+T[k-1]-Q*B[k];
}
S[0]=-P*T[0]-Q*B[0];
D2=0;
C=0;
G=0;
for(i=0;i<N;i++)
{
Q=S[j];
for(k=j-1;k>=0;k--)
Q=Q*((*X)[i]-Z)+S[k];
D2=D2+Q*Q;
C=(*Y)[i]*Q+C;
G=((*X)[i]-Z)*Q*Q+G;
}
C=C/D2;
P=G/D2;
Q=D2/D1;
D1=D2;
(*A)[j]=C*S[j];
T[j]=S[j];
for(k=j-1;k>=0;k--)
{
(*A)[k]=C*S[k]+(*A)[k];
B[k]=T[k];
T[k]=S[k];
}
}
return TRUE;
}使用示例：
曲线：y=2x^2+1 CDoubleArray dx;
CDoubleArray dy;
dx.Add(-1.0);
dx.Add(0.0);
dx.Add(1.0);
dy.Add(3.0);
dy.Add(1.0);
dy.Add(3.0);
CDoubleArray dr;
dr.Add(0.0);
dr.Add(0.0);
dr.Add(0.0);
CalculateCurveParameter(&dx, &dy, 3, 3, &dr); TRACE("%f,%f,%f", dr[0], dr[1], dr[2]);
输出为1,0,2
to ：  ciweiyu(刺猬鱼)
的确项数越多越逼近，但运算量呈平方关系上升，只是希望有个通用的方法，求的任意项数的曲线表达，但给出的坐标点数量与项数相同。
你说：给出的坐标点数量与项数相同正如前面我说的：如果是n个坐标点，用n项也就是n-1次多项式就能插值每个点了，系数通过求解一个线性方程组就可以了；这就不是逼近问题了，直接就插值了；我想到的方法：方法一：线性方程组
假设n个点为Q[n], n个未知系数P[n],要求的多项式曲线为
P(t)=P[0]+P[1]*x+P[2]*x^2+…+P[n-1]*t^(n-1)
令P(i)=Q[i],i=0,2,...,n-1,于是得到一个线性方程组，形式为
A*P=Q
其中A为n*n矩阵，求解得到P方法二：Bezier函数
假设要求的多项式曲线为
P(t) =c[0]*B_(0)(t)+c[1]*B_(1)(t)+c[2]*B_(2)(t)+…+c[n-1]*B_(n-1)(t)   t in [0,1]
就是
P(t) = sum_i(c[i]*B_(i)(t))
其中c[0]...c[n-1]为待求系数，B_(i)(t)为n阶的Bezier函数；
利用bezier函数的特点，马上得到
c[0] = Q[0],   c[n-1] = Q[n-1]
P(t)在t=0，t=1处分别求导数，马上得到c[1]，c[n-2]
...
因此，整个过程只要在端点求导数就可以，端点的j阶导数只和临近的j+1点有关系，而bezier函数端点导数都是已知的，这样可以快速求得所有c[i];在c[i]出来后，如果需要转化为方法一的形式，可以利用bezier函数和1，t，t^2,...,t^(n-1)之间转化矩阵，直接求得P[n]
多项式拟合可以转化为解N元一次方程组的问题，
这个方程组叫做法方程组
假设我们要将若干个数据拟合成3次多项式
那么法方程组有四元（N+1)怎么求法方程组你应该会吧。所以重点是解N元方程组
下面是解N元一次方程组的函数最后解出来的解Solve[0],Solve[1],Solve[2]....
就是拟合出来的多项式的常数项，一次项，二次项...系数
我以前自己写的  给你参考一下吧。。int SolveLEG(double ** CoefMatrix , int ElmtCount , double ** Solve)
{
        //
        // 解多元线性方程组函数 Solve Linear Equation Group.
        //        //
        // [in]         方程组系数矩阵  CoefMatrix[][]
        // [in]         方程组元数      ElmtCount
        // [out]        方程组的解      Solve[] *
        //
        // [return]     返回方程组是否有唯一解
         //        // 计数器
         int i,j,_i;        // 标志，用于判断秩
         bool flag,isFullRank;        // 初等变换中使用的比值u
        double u;        //
        // 一、检查矩阵的行列
         //
        isFullRank = 1;
        for (i=0;i<ElmtCount;i++)
        {
                flag = 0;
                for (j=0;j<ElmtCount;j++)
                {
                        flag = flag || (CoefMatrix[i][j]!=0);
                        if (flag == 1) break;
                }
                isFullRank = isFullRank && flag;
                if (isFullRank == 0) break;
        }
        for (i=0;i<ElmtCount;i++)
        {
                flag = 0;
                for (j=0;j<ElmtCount;j++)
                {
                        flag = flag || (CoefMatrix[j][i]!=0);
                        if (flag == 1) break;
                }
                isFullRank = isFullRank && flag;
                if (isFullRank == 0) break;
        }
        if (!isFullRank)
        {
                *Solve = NULL;
                return 0;
        }        //
        // 二、逐行进行操作
         //
        for (i=0;i<ElmtCount;i++)
        {
                //
                // 若主对角线上元素是0
                //
                if (CoefMatrix[i][i]==0)
                {
                        //
                        // 搜寻可以进行交换的行，_i 目标行标，i 需交换的行标
                           //
                        for (_i=i+1;_i<ElmtCount;_i++)
                        {
                                // 判断可交换性
                                     if (CoefMatrix[_i][i]==0) continue;                                // 逐列交换，j 列标
                                     for (j=0;j<=ElmtCount;j++)
                                {
                                        double Temp;
                                        Temp=CoefMatrix[i][j];
                                        CoefMatrix[i][j]=CoefMatrix[_i][j];
                                        CoefMatrix[_i][j]=Temp;
                                }
                        }
                }                //
                // 若主对角线上元素仍为0
                //
                if (CoefMatrix[i][i]==0)
                {
                        *Solve = NULL;
                        return 0;
                }                //
                // 初等变换
                  // 作用：逐行进行计算，_i 为循环计数器，使第i列中除了第i行以外的元素变为0。
                  //
                for (_i=0;_i<ElmtCount;_i++)
                {
                        if (_i==i) continue;                        // 计算比值
                        u = -(CoefMatrix[_i][i]/CoefMatrix[i][i]);                        // 对第_i行的逐列进行初等变换计算，j 列标
                        for (j=i;j<=ElmtCount;j++)
                                CoefMatrix[_i][j] += (CoefMatrix[i][j]*u);
                }
        }        //
        // 三、变为单位矩阵
         //
        for (i=0;i<ElmtCount;i++)
        {
                u=CoefMatrix[i][i];
                CoefMatrix[i][i] = 1;
                CoefMatrix[i][ElmtCount] /= u;
        }        //
        // 四、得到方程组的解
         //        //
        // 此时，矩阵最后一列为方程组的解         // 为方程的解开辟内存空间
         *Solve = new double[ElmtCount];        // 将最后结果存入Solve数组
         for (i=0;i<ElmtCount;i++)
                (*Solve)[i]=CoefMatrix[i][ElmtCount];        // 返回1，方程组有唯一解
         return 1;
}
http://download.csdn.net/source/2872802~~
最后还是担心你不知道怎么求法方程组求法方程组的方法：        // Power 是要拟合的次数
         // Known_x[]是已知的x序列
         // Known_y[]是已知的y序列
         // NLEGCoefMatrix[][]就是法方程组的系数矩阵了, 求到了过后就看楼上吧.
        // 注意：NLEGCoefMatrix的大小是 Power+1 * Power+2 ，不要声明成数组，用二级指针        for (i=0;i<=Power;i++)
        {
                for (j=0;j<=Power;j++)
                {
                        double Sum = 0;
                        for (k=0;k<n;k++)
                        {
                                if (i+j==0)
                                        Sum = Sum + 1;
                                else
                                        Sum = Sum + pow(Known_x[k],(double)(i+j));
                        }
                        NLEGCoefMatrix[i][j] = Sum;
                }
        }
        for (i=0;i<=Power;i++)
        {
                double Sum = 0;
                for (k=0;k<n;k++)
                {
                        if (i==0)
                                Sum = Sum + Known_y[k];
                        else
                                Sum = Sum + Known_y[k]*pow(Known_x[k],(double)i);
                }
                NLEGCoefMatrix[i][Power+1] = Sum;
        }
最后再说一句，我的答案是正确的，因为我经常用这个程序。
上面的代码有个变量没解释
n ：已知x序列和已知y序列的个数，也就是需要拟合的点的个数
http://download.csdn.net/source/2872802~~
    // 三、变为单位矩阵
         //
        for (i=0;i<ElmtCount;i++)
        {
                u=CoefMatrix[i][i];
                CoefMatrix[i][i] = 1;
                CoefMatrix[i][ElmtCount] /= u;
        }-------------顶