西安石油大学计算方法实验3--曲线拟合的最小二乘法

1、计算方法实验报告二级学院： 计算机学院 专 业： 计算机科学与技术 指导教师： 爨莹 班级学号： 201107010122 姓 名： 张文江 实验三 曲线拟合的最小二乘法1、 实验目的：在科学研究与工程技术中，常常需要从一组测量数据出发，寻找变量的函数关系的近似表达式，使得逼近函数从总体上与已知函数的偏差按某种方法度量能达到最小而又不一定过全部的点。这是工程中引入最小二曲线拟合法的出发点。 充分掌握：1最小二乘法的基本原理；2用多项式作最小二乘曲线拟合原理的基础上，通过编程实现一组实验数据的最小二乘拟合曲线。 2、实验要求: 1) 认真分析题目的条件和要求，复习相关的理论知识，选择适当的解决方

2、案和算法；2) 编写上机实验程序，作好上机前的准备工作；3) 上机调试程序，并试算各种方案，记录计算的结果（包括必要的中间结果）；4) 分析和解释计算结果；5) 按照要求书写实验报告；3、实验内容：1) 给定数据如下： x ： 0.15，0.4，0.6 ，1.01 ，1.5 ，2.2 ，2.4， 2.7， 2.9， 3.5 ，3.8 ，4.4， 4.6 ，5.1 ，6.6， 7.6； y ： 4.4964，5.1284，5.6931 ，6.2884 ，7.0989 ，7.5507 ，7.5106， 8.0756， 7.8708，8.2403 ，8.5303 ，8.7394， 8.9981 ，9

3、.1450 ，9.5070，9.9115；试作出幂函数拟合数据。 2) 已知一组数据： x ： 0，0.1，0.2 ，0.3 ，0.4 ，0.5 ，0.6，0.7，0.8，0.9 ，1 y ： -0.447，1.978，3.28 ，6.16 ，7.08 ，7.34 ，7.66，9.56，9.48，9.30 ，11.2；试用最小二乘法求多项式函数，使与此组数据相拟合。4、题目： 曲线拟合的最小二乘法 5、原理：从整体上考虑近似函数 同所给数据点 (i=0,1,m)误差 (i=0,1,m)的大小，常用的方法有以下三种：一是误差 (i=0,1,m)绝对值的最大值 ，即误差 向量 的范数；二是误差绝对

4、值的和 ，即误差向量r的1范数；三是误差平方和 的算术平方根，即误差向量r的2范数；前两种方法简单、自然，但不便于微分运算 ，后一种方法相当于考虑 2范数的平方，因此在曲线拟常采用误差平方和 来 度量误差 (i=0，1，m)的整体大小.。数据拟合的具体作法是：对给定数据 (i=0,1,，m)，在取定的函数类 中,求 ,使误差 (i=0,1,m)的平方和最小 。6、设计思想： 从几何意义上讲，就是寻求与给定点 (i=0,1,m)的距离平方和为最小的曲线。函数 称为拟合 函数或最小二乘解，求拟合函数 的方法称为曲线拟合的最小二乘法。7、对应程序：（1）幂函数程序#include<stdio.

5、h>#include<math.h>void main()double a0162,a1216,A22,Y2;int i,j,k;double x16=0.15,0.4,0.6,1.01,1.5,2.2,2.4,2.7,2.9,3.5,3.8,4.4,4.6,5.1,6.6,7.6,y16=4.4964,5.1284,5.6931,6.2884,7.0989,7.5507,7.5106,8.0756,7.8708,8.2403,8.5303,8.7394,8.9981,9.1450,9.5070,9.9115; double m0,m1,n; for(i=0;i<=15

6、;i+)a0i0=1;a0i1=log(xi);yi=log(yi);printf("输入X的值： ");printf("%f ",a0i1);printf("得到的对应的函数值：");printf("%f n ",yi);printf("n");for(i=0;i<=15;i+)for(j=0;j<=1;j+)a1ji=a0ij;/以上正确 A00=0;A01=0;A10=0;A11=0;Y0=0;Y1=0;for(i=0;i<=1;i+)for(j=0;j<=1;j+

7、) for(k=0;k<=15;k+)Aij+=a1ik*a0kj;for(i=0;i<=1;i+)for(j=0;j<=15;j+)Yi+=a1ij*yj;m0=(Y0*A11-Y1*A01)/(A00*A11-A10*A01);n=(Y0*A10-Y1*A00)/(A01*A10-A11*A00); m1=exp(m0);printf("得到的幂函数X的系数是：%fn",m1);printf("得到的幂函数X的指数是: %fn",n);（2）最小二乘法求多项式#include <stdio.h>#include <

8、conio.h>#include <math.h>#include <process.h>#define N 11/N个点#define T 3 /T次拟合#define W 1/权函数#define PRECISION 0.00001float pow_n(float a,int n)int i;if(n=0)return(1);float res=a;for(i=1;i<n;i+)res*=a;return(res);void mutiple(float aN,float bT+1,float cT+1)float res=0;int i,j,k;for

9、(i=0;i<T+1;i+)for(j=0;j<T+1;j+)res=0;for(k=0;k<N;k+)res+=aik*bkj;cij=res;void matrix_trans(float aT+1,float bN)int i,j;for(i=0;i<N;i+)for(j=0;j<T+1;j+)bji=aij;void init(float x_y2,int n)int i;printf("请输入%d个已知点：n",N);for(i=0;i<n;i+)printf("(x%d y%d):",i,i);scanf(

10、"%f %f",&x_yi0,&x_yi1);void get_A(float matrix_AT+1,float x_y2,int n)int i,j;for(i=0;i<N;i+)for(j=0;j<T+1;j+)matrix_Aij=W*pow_n(x_yi0,j);void print_array(float arrayT+1,int n)int i,j;for(i=0;i<n;i+)for(j=0;j<T+1;j+)printf("%-g ",arrayij);printf("n");

11、void convert(float arguT+2,int n)int i,j,k,p,t;float rate,temp;for(i=1;i<n;i+)for(j=i;j<n;j+)if(argui-1i-1=0)for(p=i;p<n;p+)if(argupi-1!=0)break;if(p=n)printf("方程组无解!n");exit(0);for(t=0;t<n+1;t+)temp=argui-1t;argui-1t=argupt;argupt=temp;rate=arguji-1/argui-1i-1;for(k=i-1;k<n

12、+1;k+)argujk-=argui-1k*rate;if(fabs(argujk)<=PRECISION)argujk=0;void compute(float arguT+2,int n,float root)int i,j;float temp;for(i=n-1;i>=0;i-)temp=arguin;for(j=n-1;j>i;j-)temp-=arguij*rootj;rooti=temp/arguii;void get_y(float trans_AN,float x_y2,float y,int n)int i,j;float temp;for(i=0;i&

13、lt;n;i+)temp=0;for(j=0;j<N;j+)temp+=trans_Aij*x_yj1;yi=temp;void cons_formula(float coef_AT+1,float y,float coef_formT+2)int i,j;for(i=0;i<T+1;i+)for(j=0;j<T+2;j+)if(j=T+1)coef_formij=yi;elsecoef_formij=coef_Aij;void print_root(float a,int n)int i,j;printf("%d个点的%d次拟合的多项式系数为:n",N,

14、T);for(i=0;i<n;i+)printf("a%d=%g,",i+1,ai);printf("n");printf("拟合曲线方程为:ny(x)=%g",a0);for(i=1;i<n;i+)printf(" + %g",ai);for(j=0;j<i;j+)printf("*X");printf("n");void process()float x_yN2,matrix_ANT+1,trans_AT+1N,coef_AT+1T+1,coef_for

15、muT+1T+2,yT+1,aT+1;init(x_y,N);get_A(matrix_A,x_y,N);printf("矩阵A为：n");print_array(matrix_A,N);matrix_trans(matrix_A,trans_A);mutiple(trans_A,matrix_A,coef_A);printf("正定矩阵为：n");print_array(coef_A,T+1);get_y(trans_A,x_y,y,T+1);cons_formula(coef_A,y,coef_formu);convert(coef_formu,T+1);compute(coef_f