多项式拟合

lc2236 2007-07-24 11:31:17
rt
已知若干坐标点，求出逼近曲线，即得到曲线的多项式表达式，项数不限
exam: y=c0*x0+c1*x1+c2*x2........

将如何入手？ 请做过的给予参考或提示 不胜感激
...全文
1949 12 打赏 收藏 转发到动态 举报
写回复
用AI写文章
12 条回复
切换为时间正序
请发表友善的回复…
发表回复
luo850310 2012-01-31
  • 打赏
  • 举报
 回复
// 三、变为单位矩阵
//
for (i=0;i<ElmtCount;i++)
{
u=CoefMatrix[i][i];
CoefMatrix[i][i] = 1;
CoefMatrix[i][ElmtCount] /= u;
}

-------------顶
hsejin 2011-03-29
  • 打赏
  • 举报
 回复
http://download.csdn.net/source/2872802~~
Evap 2010-12-31
  • 打赏
  • 举报
 回复
上面的代码有个变量没解释
n ：已知x序列和已知y序列的个数，也就是需要拟合的点的个数
Evap 2010-12-31
  • 打赏
  • 举报
 回复

最后还是担心你不知道怎么求法方程组

求法方程组的方法：





        // 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;

        }



最后再说一句，我的答案是正确的，因为我经常用这个程序。

hastings 2010-12-31
  • 打赏
  • 举报
 回复
http://download.csdn.net/source/2872802~~
Evap 2010-12-31
  • 打赏
  • 举报
 回复

多项式拟合可以转化为解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;

}
战争迷雾 2010-12-30
  • 打赏
  • 举报
 回复
继续搜索查找，
ciweiyu 2007-07-26
  • 打赏
  • 举报
 回复
你说：给出的坐标点数量与项数相同

正如前面我说的：如果是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]


lc2236 2007-07-25
  • 打赏
  • 举报
 回复
to ： ciweiyu(刺猬鱼)
的确项数越多越逼近，但运算量呈平方关系上升，只是希望有个通用的方法，求的任意项数的曲线表达 ，但给出的坐标点数量与项数相同。
cpio 2007-07-25
  • 打赏
  • 举报
 回复
最小二乘法曲线拟合：


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
thestral 2007-07-25
  • 打赏
  • 举报
 回复
比较好的方法是，利用正交多项式族进行拟和，比较常用的是勒让德多项式，这个在一些数值分析的书上有具体算法，比较容易历届但是步骤繁琐。对于高次的多项是拟和来说插值是不能用的，因为插值是一个高次不稳定的逼近方法。
ciweiyu 2007-07-25
  • 打赏
  • 举报
 回复
1. 已知n个坐标点，求出逼近曲线，就是说要用一个多项式曲线拟合这些离散点咯，有很多使用样条曲线进行拟合的
2. 得到曲线的多项式表达式，项数不限，这个问题就模糊了，项数越多，越逼近，如果是n个坐标点，用n项也就是n-1次多项式就能插值每个点了，系数通过求解一个线性方程组就可以了；不过，次数太高，不知意义何在
多项式拟合_多项式拟合_VBa_
利用excel vba进行编程，实现了多项式拟合
legendre_勒让德多项式拟合_多项式拟合_
使用勒让德多项式拟合函数，可调节多项式阶数
多项式拟合c++程序
对数据进行多项式拟合，可得拟合多项式的系数，法方程的系数矩阵，平方误差。
多项式拟合的C实现
多项式拟合的C实现,Linux下使用Gcc或者g++编译
C++实现最小二乘法一元回归和多项式拟合
此类的具体使用方法见：http://blog.csdn.net/czyt1988/article/details/21743595 在进行曲线拟合时用的最多的是最小二乘法，其中以一元函数（线性）和多元函数（多项式）居多，下面这个类专门用于进行多项式拟合，可以根据用户输入的阶次进行多项式拟合，算法来自于网上，和GSL的拟合算法对比过，没有问题。此类在拟合完后还能计算拟合之后的误差：SSE(剩余平方和)，SSR(回归平方和)，RMSE(均方根误差)，R-square(确定系数)。
图形处理/算法

19,468

社区成员

50,698

社区内容

发帖
与我相关
我的任务
社区描述
VC/MFC 图形处理/算法
社区管理员
  • 图形处理/算法社区
加入社区
  • 近7日
  • 近30日
  • 至今

加载中

查看更多榜单
社区公告
暂无公告

试试用AI创作助手写篇文章吧

+ 用AI写文章