分式多项式拟合VC类及应用实例

1、本文中实现了多项分式拟合类，形式如下所示（包含了多项式拟合）：y=(a0+a1x+a2xx+a3xxx.)/(1+b0x+b1xx+b2xxx.)，并且可以断项拟合，例如：a+bxx这样没有次项的方程同样能够依次输出a、b的值。通过最小二乘和迭代两种方式进行拟合。CDataFit.h文件/文件名称：CDataFit.h/文件描述：数据拟合类/作 者：dragon /完成日期：-9-26/末次修改：-9-26/本程序为数据拟合类，拟合方程形式为y=(a0+a1x+a2xx+a3xxx.)/(1+b0x+b1xx+b2xxx.)，包含了多项式拟合，并且可以断项拟合，例如：a+bxx这样没有次项的方

2、程同样能够依次输出a、b的值/拟合过程全部为C语言（VC）实现/分为最小二乘拟合和迭代拟合两种方式，分两个函数给出,两种方式在在方程有唯一解，并且迭代精度设置足够的情况下结果相同。/若不调用函数SetNumeratorDenominator设定函数形式，则默认拟合形式为y=ax/(1+bx);/#pragma onceclass CDataFitpublic:CDataFit(void);CDataFit(void);private:/确定分式形式的参数int m_nNumeratorNum;/分子项数int m_nDenominatorNum;/分母未知数n 分母项数为n+1int *m_p

3、Numerator;/指向存放分子各项的有效标记分别对应a0、a1。只能为或否则出错，表示无此项，表示有此项int *m_pDenominator;/指向存放分母各项的有效标记b0、b1。只能为或否则出错private:/解方程int m_nn;/用于存放方程中未知数的个数double* m_fDataRelative;/函数SolveResult：解方程/pRelativeCMatrix：方程系数增广矩阵例如:对于方程组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：为方程个数/n：为未知数个数/返回值：方程组有唯一解1：无穷多解2：方程组无解-1：最小二乘解/int SolveResult(double *pRelativeCMatrix,double *pResult,int m ,int n);/函数SolveResultDieDai：迭代解方程/head：方程系数增广矩阵（维数组mn+1首地址）例如:对于方程组a1*x1+a2*x2+a3*x3 = c 则head=a1,a2,a3,c/x：迭代初值，维数为n 及返回结果/m:方程个数/n：未知数个数/nn:返回迭代次数/tj:最终迭代误差维数为n/e:迭代精度/Num

5、Max:最大迭代次数默认为/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:/函数SetNumeratorDenominator：设置方程形式的函数/pNumerator：分子标记项；只能或/pDenominator：分母标记项；只能或/nNumerator

6、Num：a的总数量/nDenominatorNum：b的总数量/分式形式为：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);/函数GetDataFitRes

7、ult：获得最小二乘拟合结果/pDataX：X向量/pDataY：Y向量/n：拟合使用点数/a：拟合结果依次为分子升幂排列后的非零项系数an，分母的非零项系数bn,零项（即断项系数为）系数不输出例如y=ax/(1+bx);只输出a和b的值，分子的常数项为不再输出/返回值：成功，其他失败/int GetDataFitResult(double *pDataX,double *pDataY,int n,double *a);/函数GetDataFitIterate：获得迭代拟合结果/pDataX：X向量/pDataY：Y向量/n：拟合使用点数/a：初始输入向量，并返回拟合结果依次为a0.an,b0

8、.bn/a依次为分子升幂排列后的非零项系数an，分母的非零项系数bn,零项（即断项系数为）系数不输出例如y=ax/(1+bx);只输出a和b的值，分子的常数项为不再输出/nn:返回迭代次数/tj:最终迭代误差维数为n/e:迭代精度/NumMax:最大迭代次数默认为/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;/分子项数m_nDenominatorNum=1;/分母未知数n 分母项数为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;/分子项数m_nDenominatorNum= nDenominatorNum;/分母未知数n 分母项数为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)/获取方程失败！return -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+)/填充系数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;/函数SolveResult：解方程/pRelativeCMatrix：方程系数增广矩阵例如:对于方程组a1*x1+a2*x2+a3*x3 = c 则pRelative=a1,a2,a3,c/pResult：为方程的解，若：pRelative=a1,a2,a3,c 则pResult=x1,x2,x3;/m：为方程个数/n：为未知数个数/返回值：方程组有唯一解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; /以下是找出系数矩阵与增广矩阵的秩，比较大小判断解的情况int xishu_zhi=0,zengguang_zhi=0;/求系数矩阵秩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+)/清零，用于存放最小二乘转换后的系数曾广矩阵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+)/每次都先恢复系数矩阵memcpy(ju_zhen1,pRelativeCMatrix,sizeof(double)*nnNum);for(k=0;k<n;k+)ju_zhen1k*(n+1)+t=pRelativeCMatrixk*(n+1)+n;/再把常数列赋值给某一列/求行列式的值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;/原系数矩阵已阶梯化，求行列式简单for ( i=0;i<n;i+)pResulti=pResulti/D;/储存各未知数的解 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;/再把常数列赋值给某一列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;/原系数矩阵已阶梯化，求行列式简单for ( i=0;i<n;i+)pResul

30、ti=pResulti/D;/储存各未知数的解 delete ju_zhen1;return 0; if(xishu_zhi=zengguang_zhi && xishu_zhi<n)/第三种情况，无穷多解需求出基础解系和特解delete ju_zhen1;return 1;return 0;/函数SolveResultDieDai：迭代解方程/head：方程系数增广矩阵（维数组mn+1首地址）例如:对于方程组a1*x1+a2*x2+a3*x3 = c 则head=a1,a2,a3,c/x：迭代初值，维数为n/m:方程个数/n：未知数个数/nn:返回迭代次数/tj:最终迭

31、代误差维数为n/e:迭代精度/NumMax:最大迭代次数默认为/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