分式多項式擬合VC類及應用實例

1、本文中實現(xiàn)了多項分式擬合類，形式如下所示（包含了多項式擬合）：y=(a0+a1x+a2xx+a3xxx.)/(1+b0x+b1xx+b2xxx.)，并且可以斷項擬合，例如：a+bxx這樣沒有次項的方程同樣能夠依次輸出a、b的值。通過最小二乘和迭代兩種方式進行擬合。CDataFit.h文件/文件名稱：CDataFit.h/文件描述：數(shù)據(jù)擬合類/作 者：dragon /完成日期：-9-26/末次修改：-9-26/本程序為數(shù)據(jù)擬合類，擬合方程形式為y=(a0+a1x+a2xx+a3xxx.)/(1+b0x+b1xx+b2xxx.)，包含了多項式擬合，并且可以斷項擬合，例如：a+bxx這樣沒有次項的方

2、程同樣能夠依次輸出a、b的值/擬合過程全部為C語言（VC）實現(xiàn)/分為最小二乘擬合和迭代擬合兩種方式，分兩個函數(shù)給出,兩種方式在在方程有唯一解，并且迭代精度設置足夠的情況下結果相同。/若不調用函數(shù)SetNumeratorDenominator設定函數(shù)形式，則默認擬合形式為y=ax/(1+bx);/#pragma onceclass CDataFitpublic:CDataFit(void);CDataFit(void);private:/確定分式形式的參數(shù)int m_nNumeratorNum;/分子項數(shù)int m_nDenominatorNum;/分母未知數(shù)n 分母項數(shù)為n+1int *m_p

3、Numerator;/指向存放分子各項的有效標記分別對應a0、a1。只能為或否則出錯，表示無此項，表示有此項int *m_pDenominator;/指向存放分母各項的有效標記b0、b1。只能為或否則出錯private:/解方程int m_nn;/用于存放方程中未知數(shù)的個數(shù)double* m_fDataRelative;/函數(shù)SolveResult：解方程/pRelativeCMatrix：方程系數(shù)增廣矩陣例如:對于方程組a1*x1+a2*x2+a3*x3 = c 則pRelative=a1,a2,a3,c/pResult：為方程的解，若：pRelative=a1,a2,a3,c 則pResu

4、lt=x1,x2,x3;/m：為方程個數(shù)/n：為未知數(shù)個數(shù)/返回值：方程組有唯一解1：無窮多解2：方程組無解-1：最小二乘解/int SolveResult(double *pRelativeCMatrix,double *pResult,int m ,int n);/函數(shù)SolveResultDieDai：迭代解方程/head：方程系數(shù)增廣矩陣（維數(shù)組mn+1首地址）例如:對于方程組a1*x1+a2*x2+a3*x3 = c 則head=a1,a2,a3,c/x：迭代初值，維數(shù)為n 及返回結果/m:方程個數(shù)/n：未知數(shù)個數(shù)/nn:返回迭代次數(shù)/tj:最終迭代誤差維數(shù)為n/e:迭代精度/Num

5、Max:最大迭代次數(shù)默認為/int SolveResultDieDai(double *head,double*x,int m,int n,int &nn,double* tj,double e=0.0000000000001,int NumMax=20000);/ 迭代過程/獲得方程組int GetEquations(double *pDataX,double *pDataY,int n);public:/函數(shù)SetNumeratorDenominator：設置方程形式的函數(shù)/pNumerator：分子標記項；只能或/pDenominator：分母標記項；只能或/nNumerator

6、Num：a的總數(shù)量/nDenominatorNum：b的總數(shù)量/分式形式為：y=(pNumerator0*a0+pNumerator1*a1x+pNumerator2*a2xx+pNumerator3*a3xxx.)/(1+pDenominator0*b0x+pDenominator1*b1xx+pDenominator2*b2xxx.)，/返回值：成功，其他：失敗/int SetNumeratorDenominator(int*pNumerator,int*pDenominator,int nNumeratorNum,int nDenominatorNum);/函數(shù)GetDataFitRes

7、ult：獲得最小二乘擬合結果/pDataX：X向量/pDataY：Y向量/n：擬合使用點數(shù)/a：擬合結果依次為分子升冪排列后的非零項系數(shù)an，分母的非零項系數(shù)bn,零項（即斷項系數(shù)為）系數(shù)不輸出例如y=ax/(1+bx);只輸出a和b的值，分子的常數(shù)項為不再輸出/返回值：成功，其他失敗/int GetDataFitResult(double *pDataX,double *pDataY,int n,double *a);/函數(shù)GetDataFitIterate：獲得迭代擬合結果/pDataX：X向量/pDataY：Y向量/n：擬合使用點數(shù)/a：初始輸入向量，并返回擬合結果依次為a0.an,b0

8、.bn/a依次為分子升冪排列后的非零項系數(shù)an，分母的非零項系數(shù)bn,零項（即斷項系數(shù)為）系數(shù)不輸出例如y=ax/(1+bx);只輸出a和b的值，分子的常數(shù)項為不再輸出/nn:返回迭代次數(shù)/tj:最終迭代誤差維數(shù)為n/e:迭代精度/NumMax:最大迭代次數(shù)默認為/int GetDataFitIterate(double *pDataX,double *pDataY,int n,double *a,int &nn,double* tj,double e=0.0000000000001,int NumMax=20000);CDataFit.CPP文件#include "StdA

9、fx.h"#include "DataFit.h"#include "math.h"CDataFit:CDataFit(void)/默認值為測寬的擬合公式:y=ax/(1+bx)m_nNumeratorNum=2;/分子項數(shù)m_nDenominatorNum=1;/分母未知數(shù)n 分母項數(shù)為n+1m_pNumerator=new intm_nNumeratorNum;/指向存放分子各項的有效標記m_pDenominator=new intm_nDenominatorNum;/指向存放分母各項的有效標記m_pNumerator0=0;m_pNume

10、rator1=1;m_pDenominator0=1;/m_nn= 2;m_fDataRelative=NULL;CDataFit:CDataFit(void)if(m_pNumerator!=NULL)deletem_pNumerator;if(m_pDenominator!=NULL)deletem_pDenominator;if(m_fDataRelative!=NULL)deletem_fDataRelative;m_pNumerator= NULL;m_pDenominator= NULL;m_fDataRelative= NULL;int CDataFit:SetNumerator

11、Denominator(int*pNumerator,int*pDenominator,int nNumeratorNum,int nDenominatorNum)m_nNumeratorNum= nNumeratorNum;/分子項數(shù)m_nDenominatorNum= nDenominatorNum;/分母未知數(shù)n 分母項數(shù)為n+1if(m_pNumerator!=NULL)deletem_pNumerator;m_pNumerator=NULL;if(m_pDenominator!=NULL)deletem_pDenominator;m_pDenominator=NULL;m_pNume

12、rator= new intm_nNumeratorNum;/指向存放分子各項的有效標記m_pDenominator= new intm_nDenominatorNum;/指向存放分母各項的有效標記memcpy(m_pNumerator,pNumerator,nNumeratorNum*sizeof(int);memcpy(m_pDenominator,pDenominator,nDenominatorNum*sizeof(int);return 0;int CDataFit:GetDataFitResult(double *pDataX,double *pDataY,int n,double

13、 *a)/int i=0,j=0;if(GetEquations(pDataX,pDataY,n)!=0)/獲取方程失??！return -1;/for(i=0;i<n;i+)/TRACE(L"n");/for(j=0;j<m_nn+1;j+)/TRACE(L"%0.8ft",m_fDataRelativei*(m_nn+1)+j);/TRACE(L"n");SolveResult(m_fDataRelative,a,n,m_nn);return 0;int CDataFit:GetDataFitIterate(double

14、 *pDataX,double *pDataY,int n,double *a,int &nn,double* tj,double e,int NumMax)if(GetEquations(pDataX,pDataY,n)!=0)/獲取方程失?。eturn -1;SolveResultDieDai(m_fDataRelative,a,n,2,nn,tj,e,NumMax);return 0;int CDataFit:GetEquations(double *pDataX,double *pDataY,int n)int i=0,j=0;int nUnknownNum=0;for(i=

15、0;i<m_nNumeratorNum;i+)if(m_pNumeratori!=0)nUnknownNum+;for(i=0;i<m_nDenominatorNum;i+)if(m_pDenominatori!=0)nUnknownNum+;m_nn = nUnknownNum;if(m_fDataRelative!=NULL)deletem_fDataRelative;m_fDataRelative= NULL;m_fDataRelative = new doublen*(nUnknownNum+1);for(i=0;i<n;i+)/填充系數(shù)int lie=0;for(j

16、=0;j<m_nNumeratorNum;j+)if(m_pNumeratorj!=0)m_fDataRelativei*(nUnknownNum+1)+lie = pow(pDataXi,j);lie+;for(j=0;j<m_nDenominatorNum;j+)if(m_pDenominatorj!=0)m_fDataRelativei*(nUnknownNum+1)+lie = -pDataYi*pow(pDataXi,j+1);lie+;if(lie!=nUnknownNum)/建立方程錯誤return -1;m_fDataRelativei*(nUnknownNum+1

17、)+lie = pDataYi;/m_fDataRelativei0=pDataXi;/m_fDataRelativei1=pDataXi*pDataYi;/m_fDataRelativei2=pDataYi;return 0;/函數(shù)SolveResult：解方程/pRelativeCMatrix：方程系數(shù)增廣矩陣例如:對于方程組a1*x1+a2*x2+a3*x3 = c 則pRelative=a1,a2,a3,c/pResult：為方程的解，若：pRelative=a1,a2,a3,c 則pResult=x1,x2,x3;/m：為方程個數(shù)/n：為未知數(shù)個數(shù)/返回值：方程組有唯一解1：無窮多解

18、2：方程組無解-1：最小二乘解/int CDataFit:SolveResult(double *pRelativeCMatrix,double *pResult,int m ,int n)int i=0,j=0,t=0,k=0;int nnNum = m*(n+1);double *ju_zhen1 = new doublennNum;/求解會改變pRelative的值，因此新建一個memcpy(ju_zhen1,pRelativeCMatrix,sizeof(double)*nnNum);/把矩陣階梯化int tt=0;for(int l=0;l<n+1;l+)for(int h=t

19、t;h<m;h+)if(ju_zhen1h*(n+1)+l=0)continue; if(ju_zhen1h*(n+1)+l!=0 && h>tt)double a=0;for(int i=0;i<n+1;i+)a=-1*ju_zhen1h*(n+1)+i;ju_zhen1h*(n+1)+i=ju_zhen1tt*(n+1)+i;ju_zhen1tt*(n+1)+i=a; for(int i=tt+1;i<m;i+)if(ju_zhen1i*(n+1)+l!=0)for(int t=n;t>l;t-) ju_zhen1i*(n+1)+t=ju_zh

20、en1i*(n+1)+t-ju_zhen1i*(n+1)+l*ju_zhen1tt*(n+1)+t/ju_zhen1tt*(n+1)+l; ju_zhen1i*(n+1)+l=0;tt=tt+1; /以下是找出系數(shù)矩陣與增廣矩陣的秩，比較大小判斷解的情況int xishu_zhi=0,zengguang_zhi=0;/求系數(shù)矩陣秩int m1=0;if(m>n)m1=n; for(int i=0;i<m1;i+)for(int j=i;j<n;j+)if(ju_zhen1i*(n+1)+j!=0)xishu_zhi=xishu_zhi+1;break; /求增廣矩陣秩if(m

21、>n+1)m1=n+1; for(int i=0;i<m1;i+)for(int j=i;j<n+1;j+)if(ju_zhen1i*(n+1)+j!=0)zengguang_zhi=zengguang_zhi+1;break; if(xishu_zhi<zengguang_zhi)/第一種情況無解if(xishu_zhi=n)/超定方程求唯一最小二乘解for(i=0;i<m;i+)/清零，用于存放最小二乘轉換后的系數(shù)曾廣矩陣for(j=0;j<n+1;j+)ju_zhen1i*(n+1)+j=0.0;for(i=0;i<n;i+)/獲得最小二乘的轉換

22、矩陣for(j=0;j<n+1;j+)for(k=0;k<m;k+)ju_zhen1i*(n+1)+j += (pRelativeCMatrixk*(n+1)+j*pRelativeCMatrixk*(n+1)+i);/把矩陣階梯化int tt=0;for(int l=0;l<n+1;l+)for(int h=tt;h<m;h+)if(ju_zhen1h*(n+1)+l=0)continue; if(ju_zhen1h*(n+1)+l!=0 && h>tt)double a=0;for(int i=0;i<n+1;i+)a=-1*ju_zhe

23、n1h*(n+1)+i;ju_zhen1h*(n+1)+i=ju_zhen1tt*(n+1)+i;ju_zhen1tt*(n+1)+i=a; for(int i=tt+1;i<m;i+)if(ju_zhen1i*(n+1)+l!=0)for(int t=n;t>l;t-) ju_zhen1i*(n+1)+t=ju_zhen1i*(n+1)+t-ju_zhen1i*(n+1)+l*ju_zhen1tt*(n+1)+t/ju_zhen1tt*(n+1)+l; ju_zhen1i*(n+1)+l=0;tt=tt+1; /還原回階梯化后的值,為解最小二乘的轉換方程做準備memcpy(pRe

24、lativeCMatrix,ju_zhen1,sizeof(double)*nnNum);for(t=0;t<n;t+)/每次都先恢復系數(shù)矩陣memcpy(ju_zhen1,pRelativeCMatrix,sizeof(double)*nnNum);for(k=0;k<n;k+)ju_zhen1k*(n+1)+t=pRelativeCMatrixk*(n+1)+n;/再把常數(shù)列賦值給某一列/求行列式的值int tt=0;for(int l=0;l<n+1;l+)for(int h=tt;h<m;h+)if(ju_zhen1h*(n+1)+l=0)continue; i

25、f(ju_zhen1h*(n+1)+l!=0 && h>tt)double a=0;for(int i=0;i<n+1;i+)a=-1*ju_zhen1h*(n+1)+i;ju_zhen1h*(n+1)+i=ju_zhen1tt*(n+1)+i;ju_zhen1tt*(n+1)+i=a; for(int i=tt+1;i<m;i+)if(ju_zhen1i*(n+1)+l!=0)for(int t=n;t>l;t-) ju_zhen1i*(n+1)+t=ju_zhen1i*(n+1)+t-ju_zhen1i*(n+1)+l*ju_zhen1tt*(n+1

26、)+t/ju_zhen1tt*(n+1)+l; ju_zhen1i*(n+1)+l=0;tt=tt+1; pResultt=1;for(int i=0;i<n;i+)pResultt=pResultt*ju_zhen1i*(n+1)+i;double D=1;for(i=0;i<n;i+)D=D*pRelativeCMatrixi*(n+1)+i;/原系數(shù)矩陣已階梯化，求行列式簡單for ( i=0;i<n;i+)pResulti=pResulti/D;/儲存各未知數(shù)的解 delete ju_zhen1;return -1;elsedelete ju_zhen1;return

27、 2;if(xishu_zhi=zengguang_zhi && xishu_zhi=n)/第二種情況唯一解/利用克萊姆公式求解/還原回階梯化后的值 memcpy(pRelativeCMatrix,ju_zhen1,sizeof(double)*nnNum); for(t=0;t<n;t+)memcpy(ju_zhen1,pRelativeCMatrix,sizeof(double)*nnNum);for(k=0;k<n;k+)ju_zhen1k*(n+1)+t=pRelativeCMatrixk*(n+1)+n;/再把常數(shù)列賦值給某一列int tt=0;for(i

28、nt l=0;l<n+1;l+)for(int h=tt;h<m;h+)if(ju_zhen1h*(n+1)+l=0)continue; if(ju_zhen1h*(n+1)+l!=0 && h>tt)double a=0;for(int i=0;i<n+1;i+)a=-1*ju_zhen1h*(n+1)+i;ju_zhen1h*(n+1)+i=ju_zhen1tt*(n+1)+i;ju_zhen1tt*(n+1)+i=a; for(int i=tt+1;i<m;i+)if(ju_zhen1i*(n+1)+l!=0)for(int t=n;t>

29、;l;t-) ju_zhen1i*(n+1)+t=ju_zhen1i*(n+1)+t-ju_zhen1i*(n+1)+l*ju_zhen1tt*(n+1)+t/ju_zhen1tt*(n+1)+l; ju_zhen1i*(n+1)+l=0;tt=tt+1; pResultt=1;for(int i=0;i<n;i+)pResultt=pResultt*ju_zhen1i*(n+1)+i;double D=1;for(i=0;i<n;i+)D=D*pRelativeCMatrixi*(n+1)+i;/原系數(shù)矩陣已階梯化，求行列式簡單for ( i=0;i<n;i+)pResul

30、ti=pResulti/D;/儲存各未知數(shù)的解 delete ju_zhen1;return 0; if(xishu_zhi=zengguang_zhi && xishu_zhi<n)/第三種情況，無窮多解需求出基礎解系和特解delete ju_zhen1;return 1;return 0;/函數(shù)SolveResultDieDai：迭代解方程/head：方程系數(shù)增廣矩陣（維數(shù)組mn+1首地址）例如:對于方程組a1*x1+a2*x2+a3*x3 = c 則head=a1,a2,a3,c/x：迭代初值，維數(shù)為n/m:方程個數(shù)/n：未知數(shù)個數(shù)/nn:返回迭代次數(shù)/tj:最終迭

31、代誤差維數(shù)為n/e:迭代精度/NumMax:最大迭代次數(shù)默認為/int CDataFit:SolveResultDieDai(double *head,double*x,int m,int n,int &nn,double* tj,double e,int NumMax)/ 迭代過程/n=n+1;int i=0,k=0,j=0;double *pl=NULL,s=0.0;double *y=new doublem;/(double*)calloc(m,sizeof(double);double *t=new doublen;/(double*)calloc(n,sizeof(doubl

32、e);double *headl=new doublen;/(double*)calloc(n,sizeof(double);int b=m*(n+1);for(k=0;k<=n-1;k+)double d=0.0;for(b=0;b<=m-1;b+)d=d+headb*(n+1)+k*headb*(n+1)+k;headlk=d;for(i=0;i<=m-1;i+)s=0.0;for(j=0;j<=n-1;j+)/初始化Y0 s=s+headi*(n+1)+j*xj;yi=s;int nn = 0;for(;)/ 迭代過程for(j=0;j<=n-1;j+)tj

33、=0.0;for(k=0;k<=m-1;k+)tj=tj+headj+k*(n+1)*(head(k+1)*(n+1)-1 - yk)/headlj;for(k=0;k<=m-1;k+)yk=yk+tj*headj+k*(n+1);xj=xj+tj;nn+;for(j=0;j<=n-1;j+)/ 迭代終止的判斷if(tj>e)|(tj<-e)break;if(j=n)break;if(nn>NumMax)break;deletey;deletet;deleteheadl;memcpy(tj,t,sizeof(double)*n);if(nn>NumMa

34、x)return 1;elsereturn 0;應用程序：int i = 0,j=0;double X112 = 930.87982, 1109.5341, 1286.8435, 1465.4034, 1643.9147, 1819.9568, 2525.0876, 2700.9749, 2876.9221, 3051.9851, 3227.2148, 3400.2964;double X212 = 989.45697, 1161.9958, 1333.4760, 1506.0150, 1678.3768, 1850.5952, 2540.6509, 2715.0417, 2888.0398, 3060.8259, 3235.1245, 3407.1855;double Y12 = 0, 49.95, 99.92, 149.95, 199.90