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

1、計算方法實驗報告二級學院： 計算機學院 專 業(yè)： 計算機科學與技術 指導教師： 爨瑩 班級學號： 201107010122 姓 名： 張文江 實驗三 曲線擬合的最小二乘法1、 實驗目的：在科學研究與工程技術中，常常需要從一組測量數據出發(fā)，尋找變量的函數關系的近似表達式，使得逼近函數從總體上與已知函數的偏差按某種方法度量能達到最小而又不一定過全部的點。這是工程中引入最小二曲線擬合法的出發(fā)點。 充分掌握：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