西安石油大學計算方法實驗3--曲線擬合的最小二乘法

1、計算方法實驗報告二級學院： 計算機學院 專 業(yè)： 計算機科學與技術 指導教師： 爨瑩 班級學號： 201107010122 姓 名： 張文江 實驗三 曲線擬合的最小二乘法1、 實驗目的：在科學研究與工程技術中，常常需要從一組測量數(shù)據(jù)出發(fā)，尋找變量的函數(shù)關系的近似表達式，使得逼近函數(shù)從總體上與已知函數(shù)的偏差按某種方法度量能達到最小而又不一定過全部的點。這是工程中引入最小二曲線擬合法的出發(fā)點。 充分掌握：1最小二乘法的基本原理；2用多項式作最小二乘曲線擬合原理的基礎上，通過編程實現(xiàn)一組實驗數(shù)據(jù)的最小二乘擬合曲線。 2、實驗要求: 1) 認真分析題目的條件和要求，復習相關的理論知識，選擇適當?shù)慕鉀Q方

2、案和算法；2) 編寫上機實驗程序，作好上機前的準備工作；3) 上機調試程序，并試算各種方案，記錄計算的結果（包括必要的中間結果）；4) 分析和解釋計算結果；5) 按照要求書寫實驗報告；3、實驗內容：1) 給定數(shù)據(jù)如下： 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；試作出冪函數(shù)擬合數(shù)據(jù)。 2) 已知一組數(shù)據(jù)： 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；試用最小二乘法求多項式函數(shù)，使與此組數(shù)據(jù)相擬合。4、題目： 曲線擬合的最小二乘法 5、原理：從整體上考慮近似函數(shù) 同所給數(shù)據(jù)點 (i=0,1,m)誤差 (i=0,1,m)的大小，常用的方法有以下三種：一是誤差 (i=0,1,m)絕對值的最大值 ，即誤差 向量 的范數(shù)；二是誤差絕對

4、值的和 ，即誤差向量r的1范數(shù)；三是誤差平方和 的算術平方根，即誤差向量r的2范數(shù)；前兩種方法簡單、自然，但不便于微分運算 ，后一種方法相當于考慮 2范數(shù)的平方，因此在曲線擬常采用誤差平方和 來 度量誤差 (i=0，1，m)的整體大小.。數(shù)據(jù)擬合的具體作法是：對給定數(shù)據(jù) (i=0,1,，m)，在取定的函數(shù)類 中,求 ,使誤差 (i=0,1,m)的平方和最小 。6、設計思想： 從幾何意義上講，就是尋求與給定點 (i=0,1,m)的距離平方和為最小的曲線。函數(shù) 稱為擬合 函數(shù)或最小二乘解，求擬合函數(shù) 的方法稱為曲線擬合的最小二乘法。7、對應程序：（1）冪函數(shù)程序#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("得到的對應的函數(shù)值：");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("得到的冪函數(shù)X的系數(shù)是：%fn",m1);printf("得到的冪函數(shù)X的指數(shù)是: %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/權函數(shù)#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次擬合的多項式系數(shù)為: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