最小二乘法多项式拟合c语言

最小二乘法多项式拟合

对于给定的数据点Niyxii1),,(，可用下面的n阶多项式进行拟合，即

为了使拟合出的近似曲线能尽量反映所给数据的变化趋势，要求在所有数据点上的残差

都较小。为达到上述目标，可以令上述偏差的平方和最小，即

称这种方法为最小二乘原则，利用这一原则确定拟合多项式)(xf的方法即为最小二乘法多项式拟合。

确定上述多项式的过程也就是确定)(xf中的系数nkak0,的过程，根据最小二乘原则，则偏差平方和应该是这些系数的函数，即

为使上式取值最小，则其关于nkak0,的一阶导数应该为零，即有

将上面各等式写成方程组的形式可有

写成矩阵形式有

上述方程组可以通过克莱姆法则来计算，从而解出各系数nkak0,得到拟合方程。

考虑到一般情况提高拟合多项式的阶数并不能提高拟合精度，所以常用的多项拟合阶数为一阶和二阶，即线性拟合和二次拟合。两者的计算公式如下：

关于线性拟合，除上面按克莱姆法则来计算外，还可以有另一思路，下面对此进行说明。由于是线性拟合，最后得到的是一条直线，因此，直线可以由斜率和截距两个参数来确定，因此，求出这两个参数即可。首先对克莱姆法的求解结果进行展开可以得到

下面考虑先计算斜率再计算截距的方法，从下图可见，斜率计算与坐标系的位置无关，所以可以将坐标原点平移到样本的ix和iy坐标的均值所在点上

图中

则在新的坐标系),(yx下斜率的计算公式与前面1a的计算公式相同，将其中的坐标),(yx换成),(yx即可得到下面的计算公式

由样本在新坐标系下的坐标ix和iy的均值为零，或者由下面推导可知

则斜率的计算公式可以简化为

还原为原坐标有

下面推导截距的计算公式 xyxy),(yx这样可以得到两组计算公式，分别如下

或

/*最小二乘法的曲线拟合*/

#include

#include

#include

#definemax100

voidmain()

{

inti,j,k,m,N,mi;

floatmx,temp;

floatX[max][max],Y[max],x[max],y[max],a[max];

FILE*fp1;

if((fp1=fopen("in1.txt","r"))==NULL)/*输入拟合曲线的次数m以及已知的数据

组数N*/{

printf("Can'topenthisfile!\n");

exit(0);

}

for(i=0;i<2;i++)

fscanf(fp1,"%d%d",&m,&N);

fclose(fp1);

FILE*fp2;

if((fp2=fopen("in2.txt","r"))==NULL)/*输入已知的N组数据坐标*/

{

printf("Can'topenthisfile!\n");

exit(0);

}

for(i=0;i

fscanf(fp2,"%f%f",&x[i],&y[i]);

fclose(fp2);

for(i=0;i<=m;i++)/*由给定的点得系数矩阵*/

for(j=i;j<=m;j++)

{

temp=0;

for(k=0;k

temp=temp+pow(x[k],(i+j));

X[i][j]=temp;

X[j][i]=X[i][j];

}

for(i=0;i<=m;i++)/*由给定的点得右端矩阵列向量*/

{temp=0;

for(k=0;k

temp=temp+y[k]*pow(x[k],i);

Y[i]=temp;

}

for(j=0;j

项式的系数*/{

for(i=j+1,mi=j,mx=fabs(X[j][j]);i<=m;i++)

if(fabs(X[i][j])>mx)

{

mi=i;

mx=fabs(X[i][j]);

}

if(j

{

temp=Y[j];

Y[j]=Y[mi];

Y[mi]=temp;

for(k=j;k<=m;k++)

#include

#include

//////////////////////////////////////////////////////////////////////////////////////////

//矩阵结构体

struct Matrix

{

int m,n;//m为行数,n为列数

double **pm;//指向矩阵二维数组的指针

};

//初始化矩阵mt，并置矩阵的行为m，列为n

void InitMatrix(struct Matrix *mt,int m,int n)

{

int i;

//指定矩阵的行和列

mt->m=m;

mt->n=n;

//为矩阵分配内存

mt->pm=new double *[m];

for(i=0;i

{

mt->pm[i]=new double[n];

}

}

//用0初始化矩阵mt，并置矩阵的行为m，列为n

void InitMatrix0(struct Matrix *mt,int m,int n)

{

int i,j;

InitMatrix(mt,m,n);

for(i=0;ipm[i][j]=0.0;

}

//销毁矩阵mt

void DestroyMatrix(struct Matrix *mt)

{

int i;

//释放矩阵内存

for(i=0;im;i++)

{

delete []mt->pm[i];

}

delete []mt->pm;

}

//矩阵相乘mt3=mt1*mt2

//成功返回1，失败返回0

int MatrixMul(Matrix *mt1,Matrix *mt2,Matrix *mt3)

{

int i,j,k;

double s;

//检查行列号是否匹配

if(mt1->n!=mt2->m||mt1->m!=mt3->m||mt2->n!=mt3->n) return 0;

连续系统仿真实验报告

实验数据拟合建模

姓名：

专业：

学号：

时间：2013年5月1日

实验单元二 实验数据拟合建模

一、实验目的

1、用C语言实现最小二乘的多项式拟合和LU分解法；

2、熟练掌握最小二乘拟合和LU分解法的基本原理。

3、体会用计算机编程解决计算问题的方法。

二、需求说明

(一)、需求阐述

本次实验是要求根据已知的自变量和函数值，通过多项式拟合来分别计算2、3、4阶拟合多项式，并根据拟合结果分别计算出待求点的函数值。其中解拟合系数方程组时采用LU分解的方法计算拟合多项式的系数。

（二）、实验公式

m次拟合函数公式为：

01x=++mmaaxax

计算系数ia的方程组为：

{ 00110++...+=mmsasasat1021+11++...+=mmsasasat⋮0+112++...+=mmmmsasasat ○1

其中 ++=0=nkjjkiisx,

=0=njiiityx

所以，在编程计算时，先计算出方程组○1，再用LU分解法计算求出ia的值，即可得到拟合多项式。LU分解法的公式为：

[ 1𝑙21𝑙31 ⋮𝑙𝑛101𝑙32⋮𝑙𝑛2001⋮⋯⋯⋯⋱⋱⋯0

0⋮01]

[ 𝑢 110000 𝑢 12 𝑢 22000 𝑢 13 𝑢 23 𝑢 33⋱⋯⋯⋯⋱⋱0 𝑢 1𝑛

𝑢 2𝑛⋮ 𝑢 𝑛−1𝑛 𝑢 𝑛𝑛]

=[ a11a21a31⋮an1a12a22a32⋮an2a13a23a33⋮an3⋯⋯⋯⋱⋯a1n

a2na3n⋮ann]

其中L矩阵和U矩阵的计算公式如下：

第一步，当k=1,有： { 11= (=2,3,,n)jjuaj1111111== (=2,3,,n)iiiaaliua