最小二乘拟合C语言程序

输入文件zuixiaoercheng_input.txt中内容：

5 3

-2 -1 0 1 2

-0.1 0.1 0.4 0.9 1.6

代码部分：

#include

#include

#include

#include

#define p(k) (p+(k)*m)

#define pkx(k) (pkx+(k)*n)

double *x,*f,*p,*a,*b,*d,*pkx,*sx,*xishu;

int m,n;

void zhengjiaonihe();

void size(int *m,int *n);

void input(double x[],double f[],int m,int n);

void output(int n,double pkx[],double xishu[],double R);

double NeiJi(double *F,double *G,double *W);

void GetPk(int k);

void GetXiShu(int k);

double Wucha(int k);

void main()

{

zhengjiaonihe();

system("pause");

}

void zhengjiaonihe()

{

int k,i,j;

double t,R;

printf("\n本程序是用正交函数作最小二乘拟合\n");

size(&m,&n);

x=(double *)malloc(sizeof(double)*m);

f=(double *)malloc(sizeof(double)*m);

p=(double *)malloc(sizeof(double)*m*n);

a=(double *)malloc(sizeof(double)*n);

b=(double *)malloc(sizeof(double)*n);

d=(double *)malloc(sizeof(double)*n);

pkx=(double *)malloc(sizeof(double)*n*n);

sx=(double *)malloc(sizeof(double)*m);

xishu=(double *)malloc(sizeof(double)*n);

input(x,f,m,n);

for (i=0;i

{

for (j=0;j

pkx(i)[j]=0; // pkx 初始化

}

for (k=0;k

{

p(0)[k]=1; // P0(x)=1;

}

pkx(0)[0]=1;

d[0]=NeiJi(f,p(0),p(0))/NeiJi(p(0),p(0),p(0));

GetXiShu(1);

R=Wucha(1);

a[1]=NeiJi(x,p(0),p(0))/NeiJi(p(0),p(0),p(0)); //α1=(xP1,P1)/(P1,P1)

pkx(1)[0]=-a[1];pkx(1)[1]=1;

for (k=0;k

p(1)[k]=x[k]-a[1]; // P1(x)=(x-α1)P0(x) P0(x)=1

d[1]=NeiJi(f,p(1),p(0))/NeiJi(p(1),p(1),p(0));

GetXiShu(2);

R=Wucha(2);

for (k=1;k

{

t=NeiJi(p(k),p(k),p(0));

a[k+1]=NeiJi(p(k),p(k),x)/t; // αk+1=(xPk,Pk)/(Pk,Pk)

b[k]=t/NeiJi(p(k-1),p(k-1),p(0)); //βk=(Pk,Pk)/(Pk-1,Pk-1)

GetPk(k+1); // 计算Pk(x)在各节点的函数值

d[k+1]=NeiJi(f,p(k+1),p(0))/NeiJi(p(k+1),p(k+1),p(0));

// ak=(f,Pk)/(Pk,Pk)

pkx(k+1)[0]=-a[k+1]*pkx(k)[0]-b[k]*pkx(k-1)[0];

for (i=1;i

{ // Pk+1=(x-αk+1)Pk(x)-βkPk-1(x)

pkx(k+1)[i]=pkx(k)[i-1]-a[k+1]*pkx(k)[i]-b[k]*pkx(k-1)[i];

}

GetXiShu(k+1); // 求拟合多项式的系数

R=Wucha(k+1); //计算平方误差

}

GetXiShu(n); // 求拟合多项式的系数

R=Wucha(n); //计算平方误差

output(n,pkx,xishu,R);

}

double NeiJi(double *F,double *G,double *W) //内积函数

{

double s=0;

int i;

for (i=0;i

s+=F[i]*G[i]*W[i];

return s;

} // NeiJi end

void GetPk(int k) // 求正交多项式Pk(x)在节点 x(k) 的值

{

int i;

for (i=0;i

p

(k)[i]=(x[i]-a[k])*p(k-1)[i]-b[k-1]*p(k-2)[i];

}

void GetXiShu(int k) // 求拟合多项式的系数

{

int i,j;

double t;

for (i=0;i

{

t=0;

for (j=0;j

{

t+=d[j]*pkx(j)[i];

}

xishu[i]=t;

}

} // GetXiShu end

double Wucha(int k) // 计算平方误差

{

double r1;

int i,j;

for (i=0;i

{

sx[i]=0;

for (j=0;j

sx[i]=sx[i]+xishu[j]*pow(x[i],j); // 计算拟合函数在各节点的值

}

r1=fabs(NeiJi(f,f,p(0))-NeiJi(f,sx,p(0))); // 计算平方误差

return r1;

} // Wucha end

void input(double x[],double f[],int m,int n)

{ FILE *fp;

char *filename="zuixiaoercheng_input.txt";

int j;

float tmp;

if((fp=fopen(filename,"r"))==NULL)

{printf("can not open file.\n");

system("pause");

}

fscanf(fp,"%d",&m);

fscanf(fp,"%d",&n);

printf("结点个数m=%d,次数n=%d\n",m,n);

printf("x=\n");

for(j=0;j

{fscanf(fp,"%f",&tmp);

x[j]=tmp; /* 读取x */

printf("%f ",tmp);

}

printf("\n");

printf("f=\n");

for(j=0;j

{fscanf(fp,"%f",&tmp);

f[j]=tmp; /* 读取y */

printf("%f ",tmp);

}

fclose(fp);

}

void output(int n,double pkx[],double xishu[],double R)

{ FILE *fp;

char *filename="zuixiaoercheng_output.txt";

int k,i;

if((fp=fopen(filename,"w"))==NULL)

{printf("can not open file.\n");

system("pause");

}

printf("\n正交多项式Pk(x)为: \n");

for (k=0;k

{

{printf("P%d(x)=%lG",k,pkx(k)[0]);

fprintf(fp,"P%d(x)=%lG",k,pkx(k)[0]);

for (i=1;i

if(pkx(k)[i]>=0)

{printf("+%lGx^%d",pkx(k)[i],i);

fprintf(fp,"+%lGx^%d",pkx(k)[i],i);

}

else

{printf("%lGx^%d\n",pkx(k)[i],i);

fprintf(fp,"%lGx^%d\n",pkx(k)[i],i);

}

}

printf("\n");

fprintf(fp,"\n");

}

printf("\n平方误差为: %lf\n",R);

printf("\n拟合多项式为: \n");

printf("F(x)=%lG",xishu[0]);

fprintf(fp,"\n平方误差为: %lf\n",R);

fprintf(fp,"\n拟合多项式为: \n");

fprintf(fp,"F(x)=%lG",xishu[0]);

for (k=1;k

{

if(xishu[k]>=0)

{printf("+%lGx^%d",xishu[k],k);

fprintf(fp,"+%lGx^%d",xishu[k],k);

}

else

{printf("%lGx^%d",xishu[k],k);

fprintf(fp,"%lGx^%d",xishu[k],k);

}

}

printf("\n");

fclose(fp);

}

void size(int *m,int *n)

{FILE *fp;

char *filename="zuixiaoercheng_input.txt";

if((fp=fopen(filename,"r"))==NULL)

{printf("can not open file.\n");

system("pause");

}

fscanf(fp,"%d",m);

fscanf(fp,"%

d",n);

fclose(fp);

}

输入文件zuixiaoercheng_input.txt中内容：

5 3

-2 -1 0 1 2

-0.1 0.1 0.4 0.9 1.6

代码部分：

#include

#include

#include

#include

#define p(k) (p+(k)*m)

#define pkx(k) (pkx+(k)*n)

double *x,*f,*p,*a,*b,*d,*pkx,*sx,*xishu;

int m,n;

void zhengjiaonihe();

void size(int *m,int *n);

void input(double x[],double f[],int m,int n);

void output(int n,double pkx[],double xishu[],double R);

double NeiJi(double *F,double *G,double *W);

void GetPk(int k);

void GetXiShu(int k);

double Wucha(int k);

void main()

{

zhengjiaonihe();

system("pause");

}

void zhengjiaonihe()

{

int k,i,j;

double t,R;

printf("\n本程序是用正交函数作最小二乘拟合\n");

size(&m,&n);

x=(double *)malloc(sizeof(double)*m);

f=(double *)malloc(sizeof(double)*m);

p=(double *)malloc(sizeof(double)*m*n);

a=(double *)malloc(sizeof(double)*n);

b=(double *)malloc(sizeof(double)*n);

d=(double *)malloc(sizeof(double)*n);

pkx=(double *)malloc(sizeof(double)*n*n);

sx=(double *)malloc(sizeof(double)*m);

xishu=(double *)malloc(sizeof(double)*n);

input(x,f,m,n);

for (i=0;i

{

for (j=0;j

pkx(i)[j]=0; // pkx 初始化

}

for (k=0;k

{

p(0)[k]=1; // P0(x)=1;

}

pkx(0)[0]=1;

d[0]=NeiJi(f,p(0),p(0))/NeiJi(p(0),p(0),p(0));

GetXiShu(1);

R=Wucha(1);

a[1]=NeiJi(x,p(0),p(0))/NeiJi(p(0),p(0),p(0)); //α1=(xP1,P1)/(P1,P1)

pkx(1)[0]=-a[1];pkx(1)[1]=1;

for (k=0;k

p(1)[k]=x[k]-a[1]; // P1(x)=(x-α1)P0(x) P0(x)=1

d[1]=NeiJi(f,p(1),p(0))/NeiJi(p(1),p(1),p(0));

GetXiShu(2);

R=Wucha(2);

for (k=1;k

{

t=NeiJi(p(k),p(k),p(0));

a[k+1]=NeiJi(p(k),p(k),x)/t; // αk+1=(xPk,Pk)/(Pk,Pk)

b[k]=t/NeiJi(p(k-1),p(k-1),p(0)); //βk=(Pk,Pk)/(Pk-1,Pk-1)

GetPk(k+1); // 计算Pk(x)在各节点的函数值

d[k+1]=NeiJi(f,p(k+1),p(0))/NeiJi(p(k+1),p(k+1),p(0));

// ak=(f,Pk)/(Pk,Pk)

pkx(k+1)[0]=-a[k+1]*pkx(k)[0]-b[k]*pkx(k-1)[0];

for (i=1;i

{ // Pk+1=(x-αk+1)Pk(x)-βkPk-1(x)

pkx(k+1)[i]=pkx(k)[i-1]-a[k+1]*pkx(k)[i]-b[k]*pkx(k-1)[i];

}

GetXiShu(k+1); // 求拟合多项式的系数

R=Wucha(k+1); //计算平方误差

}

GetXiShu(n); // 求拟合多项式的系数

R=Wucha(n); //计算平方误差

output(n,pkx,xishu,R);

}

double NeiJi(double *F,double *G,double *W) //内积函数

{

double s=0;

int i;

for (i=0;i

s+=F[i]*G[i]*W[i];

return s;

} // NeiJi end

void GetPk(int k) // 求正交多项式Pk(x)在节点 x(k) 的值

{

int i;

for (i=0;i

p

(k)[i]=(x[i]-a[k])*p(k-1)[i]-b[k-1]*p(k-2)[i];

}

void GetXiShu(int k) // 求拟合多项式的系数

{

int i,j;

double t;

for (i=0;i

{

t=0;

for (j=0;j

{

t+=d[j]*pkx(j)[i];

}

xishu[i]=t;

}

} // GetXiShu end

double Wucha(int k) // 计算平方误差

{

double r1;

int i,j;

for (i=0;i

{

sx[i]=0;

for (j=0;j

sx[i]=sx[i]+xishu[j]*pow(x[i],j); // 计算拟合函数在各节点的值

}

r1=fabs(NeiJi(f,f,p(0))-NeiJi(f,sx,p(0))); // 计算平方误差

return r1;

} // Wucha end

void input(double x[],double f[],int m,int n)

{ FILE *fp;

char *filename="zuixiaoercheng_input.txt";

int j;

float tmp;

if((fp=fopen(filename,"r"))==NULL)

{printf("can not open file.\n");

system("pause");

}

fscanf(fp,"%d",&m);

fscanf(fp,"%d",&n);

printf("结点个数m=%d,次数n=%d\n",m,n);

printf("x=\n");

for(j=0;j

{fscanf(fp,"%f",&tmp);

x[j]=tmp; /* 读取x */

printf("%f ",tmp);

}

printf("\n");

printf("f=\n");

for(j=0;j

{fscanf(fp,"%f",&tmp);

f[j]=tmp; /* 读取y */

printf("%f ",tmp);

}

fclose(fp);

}

void output(int n,double pkx[],double xishu[],double R)

{ FILE *fp;

char *filename="zuixiaoercheng_output.txt";

int k,i;

if((fp=fopen(filename,"w"))==NULL)

{printf("can not open file.\n");

system("pause");

}

printf("\n正交多项式Pk(x)为: \n");

for (k=0;k

{

{printf("P%d(x)=%lG",k,pkx(k)[0]);

fprintf(fp,"P%d(x)=%lG",k,pkx(k)[0]);

for (i=1;i

if(pkx(k)[i]>=0)

{printf("+%lGx^%d",pkx(k)[i],i);

fprintf(fp,"+%lGx^%d",pkx(k)[i],i);

}

else

{printf("%lGx^%d\n",pkx(k)[i],i);

fprintf(fp,"%lGx^%d\n",pkx(k)[i],i);

}

}

printf("\n");

fprintf(fp,"\n");

}

printf("\n平方误差为: %lf\n",R);

printf("\n拟合多项式为: \n");

printf("F(x)=%lG",xishu[0]);

fprintf(fp,"\n平方误差为: %lf\n",R);

fprintf(fp,"\n拟合多项式为: \n");

fprintf(fp,"F(x)=%lG",xishu[0]);

for (k=1;k

{

if(xishu[k]>=0)

{printf("+%lGx^%d",xishu[k],k);

fprintf(fp,"+%lGx^%d",xishu[k],k);

}

else

{printf("%lGx^%d",xishu[k],k);

fprintf(fp,"%lGx^%d",xishu[k],k);

}

}

printf("\n");

fclose(fp);

}

void size(int *m,int *n)

{FILE *fp;

char *filename="zuixiaoercheng_input.txt";

if((fp=fopen(filename,"r"))==NULL)

{printf("can not open file.\n");

system("pause");

}

fscanf(fp,"%d",m);

fscanf(fp,"%

d",n);

fclose(fp);

}


相关文章

  • 神经网络在数据拟合方面的应用

    • 神经网络在数据拟合方面的应用 摘要 本文将讲述人工神经网络及其数据拟合中的应用.人工神经网络是从信息处理角度对人脑神经元网络进行抽象,建立某种简单模型,按不同的连接方式组成不同的网络.它在模式识别.智能机器人.自动控制.预测估计.生物.医学 ...查看


  • 最小二乘法

    • 最小二乘法 设(x 1, y 1 ), (x 2, y 2), -, (x n, y n)是直角平面坐标系下给出 的一组数据,若x 1<x 2<-<x n,我们也可以把这组数据看作 是一个离散的函数.根据观察,如果这组数据 ...查看


  • 传感器综合实验仿真报告

    • 综合实验报告 ( 2015 -- 2016年度第一学期) 名 称: 传感器原理与应用 题 目: 院 系: 控制与计算机工程 班 级: 测控1303 学 号: 学生姓名:指导教师: 设计周数: 一周 成 绩: 日期:2016 年1月15 日 ...查看


  • 最小二乘法分段直线拟合_田垅

    • 第39卷 第6A期2012年6月计算机科学 ComutercienceSVol．39No．6A June2012最小二乘法分段直线拟合 田 垅 刘宗田 ()上海大学计算机工程与科学学院 上海200072 摘 要 曲线拟合是图像分析中非常重要 ...查看


  • Matlab数据拟合程序

    • 课程设计名称: 设计二:数据拟合 指导教师: 张莉 课程设计时数: 6 课程设计设备:安装了Matlab .C ++软件的计算机 课程设计日期: 实验地点: 第五教学楼北902 课程设计目的: 1. 了解最小二乘拟合的原理,掌握用MA TL ...查看


  • 线性拟合法

    • 摘 要 摘 要 插值法和曲线拟合是两种来源于实际,同时又广泛应用于实际的重要的数值 计算方法.随着计算机技术的不断发展以及人类计算机水平的逐步提高,他们在国民经济和科学研究中占据了越来越重要的地位.插值法与曲线拟合结合计算机技术例如MATL ...查看


  • 击实试验数据的数值分析方法

    • 击实试验数据的数值分析方法 俞文生 蒲 华2 (1 江西交通工程咨询监理中心 南昌 330008) (2 江西省交通质量监督站 南昌 330008) 1 摘 要:对于公路填方路基工程,土的最大干密度ρdm 和最佳含水量ω0是路基施工质量控制 ...查看


  • 太阳影子定位

    • 太阳影子定位 摘 要 本文考虑了影响直杆影长的五个因素--直杆所在地经度.纬度.日期(太阳积日).时间.直杆高度五个参数,引入太阳高度角.时角.太阳方位角.赤纬四个辅助变量,利用三角函数关系分析了五个参数与影长变化的关系,建立出影长变化的数 ...查看


  • 光电效应测普朗克常数实验及数据处理方法ly2010

    • 光电效应测普朗克常数实验及数据处理方法 (东南大学 仪器科学与工程学院, 南京 210096) 摘 要: 简述了用作图法处理光电效应实验数据的优点和不足,通过实例说明,采用最小二乘拟合算法处理实验数据更具客观合理性. 关键词: 光电效应: ...查看


热门内容