函数逼近与拟合法

第五讲函数迫近与拟正当第1页内容提要引言函数迫近傅里叶迫近最小二乘法拟合最小二乘法多元线性拟合非线性拟合MATLAB拟合函数小结2025/4/132第2页例：考查某种纤维强度与其拉伸倍数关系,下表是实际测定24个纤维样品强度与对应拉伸倍数统计:纤维强度随拉伸倍数增加而增加。1、引言2025/4/133第3页24个点大致分布在一条直线附近。故可认为强度y与拉伸倍数x主要关系应为线性关系：必须找到一个度量标准来衡量什么曲线最靠近全部数据点。2025/4/134第4页2025/4/135第5页在一个包含有很多数据点区间内结构插值函数，必定使用高次多项式。而高次插值多项式是不稳定。因为数据本身存在误差，利用插值方法得到插值多项式必定保留了全部测量误差，造成插值函数与物理规律差异较大。

试验数据拟合能够克服插值方法在处理这类问题中存在缺点。

对这么数据采取上一讲介绍插值方法近似求描述物理规律解析函数，必定存在以下缺点：2025/4/136第6页试验数据拟合基本思想：

使近似函数尽可能靠近数据点，而不要求近似函数一定经过全部数据点。

试验数据拟合能够在一定精度内找出反应物理量间客观函数关系解析式。假如试验数据存在误差，这种做法能够部分抵消原来数据中测量误差，从而使所得到拟合函数更加好地反应物理规律。2025/4/137第7页利用拟合能够处理两类物理问题：物理规律已知，但描述物理规律解析式中一些系数未知，能够利用试验方法取得了物理量之间关系，经过拟合方法，求出这些系数近似值。物理规律未知，利用试验方法取得了物理量之间关系，经过拟合方法，得到一个近似解析式，用于描述物理规律。拟合函数尽可能靠近数据点怎样实现？2025/4/138第8页2、函数迫近

在区间[a,b]上已知一连续函数f(x)，假如该函数表示式太过于复杂不利于进行计算机运算，就会利用一个简单函数去近似f(x)，这就是函数迫近问题。假如f(x)表示式未知，只知道描述f(x)一条曲线，这就是曲线拟合问题。和插值问题不一样，迫近和拟合并不要求迫近函数在已知点上值一定等于原函数函数值，而是按照某种标准使得二者差值到达最小。2025/4/139第9页迫近方法：Chebyshev(切比雪夫)迫近：连续函数，多项式。F=Chebyshev(y,k,x0)Legendre(勒让德)迫近：多项式。F=Legendre(y,k,x0)Pade(帕德)迫近：有理分式。F=Pade(y,k,x0)傅里叶迫近：周期函数，三角多项式。连续周期函数，[A0,A,B]=FZZ(func,T,n)离散周期函数，c=DFF(f,N)2025/4/1310第10页Chebyshev(切比雪夫)迫近当一个连续函数定义在区间[-1，1]上时，能够展开成为切比雪夫级数。2025/4/1311实际应用中，依据精度要求来截取有限项数第11页functionf=Chebyshev(y,k,x0)%用切比雪夫多项式迫近已知函数%已知函数：y%迫近已知函数所需项数：k%迫近点x坐标：x0%求得切比雪夫迫近多项式或在x0处迫近值：fsymst;T(1:k+1)=t;T(1)=1;T(2)=t;c(1:k+1)=0.0;c(1)=int(subs(y,findsym(sym(y)),sym('t'))*T(1)/sqrt(1-t^2),t,-1,1)/pi;c(2)=2*int(subs(y,findsym(sym(y)),sym('t'))*T(2)/sqrt(1-t^2),t,-1,1)/pi;f=c(1)+c(2)*t;……fori=3:k+1T(i)=2*t*T(i-1)-T(i-2);c(i)=2*int(subs(y,findsym(sym(y)),sym('t'))*T(i)/sqrt(1-t^2),t,-1,1)/2;f=f+c(i)*T(i);f=vpa(f,6);if(i==k+1)if(nargin==3)f=subs(f,'t',x0);elsef=vpa(f,6);endendend2025/4/1312第12页2025/4/1313例：用切比雪夫公式（取6项）迫近函数1/（2-x），并求当x=0.5时函数值函数准确值为1/（2-0.5）=0.6667，可见迫近结果比较靠近第13页离散周期函数傅里叶迫近functionc=DFF(f,N)%用傅里叶级数迫近已知离散周期函数%离散数据点：f%展开项数：N%离散傅里叶迫近系数：c

c(1:N)=0;for(m=1:N)for(n=1:N)c(m)=c(m)+f(n)*exp(-i*m*n*2*pi/N);endc(m)=c(m)/N;end2025/4/1314第14页例N123456Y0.84150.90930.1411-0.7568-0.9589-0.2794>>y=[0.84150.90930.1411-0.7568-0.9589-0.2794];>>c=DFF(y,6)c=Columns1through4-0.0926-0.5003i-0.0260-0.0194i-0.0251+0.0000i-0.0260+0.0194iColumns5through6-0.0926+0.5003i-0.0172-0.0000i2025/4/1315第15页3.1最小二乘法

首先，从一个简单例子来讨论一元线性拟合与最小二乘法问题。为了含有普通性，把上式改写为：

经过试验测量，求金属铜电阻温度系数α，金属电阻与温度关系以下：3、最小二乘法拟合2025/4/1316第16页

经过试验测得金属铜温度x与电阻y数据以下：xi(℃)Yi(Ω)xi(℃)Yi(Ω)xi(℃)Yi(Ω)04.38705.581406.74104.56805.741506.94204.70905.961607.12304.861006.061707.28405.081106.261807.42505.241206.441907.60605.401306.582007.782025/4/1317第17页

设一元线性拟合函数为：将试验数据代入拟合函数，得到方程组21个线性方程矛盾方程组2025/4/1318第18页

因为以上矛盾方程组不能确定一组唯一A0和A1，也就是说，由方程组可求得A0和A1多组解，那么终究哪一组解最靠近客观真实值呢?

按照拟合思想，应该使在每一个测量点拟合函数函数值尽可能靠近测量值，这么拟合函数才是满足要求，即：定义偏差：2025/4/1319第19页

按照拟合思想，必须在每一个测量点偏差都很小，怎样到达这一要求？

不过因为偏差有正有负，求和时可能相互抵消，这并不能确保在每一个测量点偏差都很小。方法一：偏差之和最小

尽管这种方法能够确保在每一个测量点偏差都很小，但这种方法数学处理比较困难。方法二：偏差绝对值之和最小2025/4/1320第20页

这种方法既能够确保在每一个测量点偏差都很小，又方便数学处理，所以这种方法是可行。方法三：偏差平方和最小-----最小二乘法2025/4/1321第21页残差向量各分量平方和记为：最小二乘法：以残差平方和最小问题解来确定拟合函数方法。令－－在回归分析中称为残差(i=1,2,…m)残差向量：2025/4/1322第22页由多元函数求极值必要条件，有可得即2025/4/1323第23页上式为由n+1个方程组成方程组，称正规方程组。由得即2025/4/1324第24页引入记号则由内积概念可知显然内积满足交换律正规方程组便可化为2025/4/1325第25页将其表示成矩阵形式：其系数矩阵为对称阵。所以正规方程组系数矩阵非奇异,即依据Crame法则,正规方程组有唯一解，称其为最小二乘解。2025/4/1326第26页2025/4/1327function[A,b,p]=Least_square(wfun,phifun,x,y)%最小二乘拟合%输入参数：%---wfun：权系数%---phifun：拟合基函数%---x,y：已知数据x,y坐标%---n：数据拟合次数，默认值为1%输出参数：%---A：法方程组系数矩阵%---b：法方程组右端向量%---p：最小二乘拟合系数phifun=phifun(x);n=size(phifun,1);A=zeros(n);fori=1:nforj=1:nA(i,j)=sum(wfun(i).*phifun(i,:).*phifun(j,:));endb(i)=sum(wfun(i).*phifun(i,:).*y);endb=b';a=A\b;p=a';MATLAB实现第27页2025/4/1328x=0:0.5:3;%x轴数据y=[00.47940.84150.98150.91260.59850.1645];%y轴数据wfun=ones(1,6);%权系数phifun=@(x)[ones(size(x));x;x.^2;cos(x);exp(x);sin(x)];%拟合基函数[A,b,p]=Least_square(wfun,phifun,x,y)%最小二乘拟合求解symsxdigits(4)%设定精度Phifun=[1;x;x.^2;cos(x);exp(x);sin(x)];y=vpa(p*Phifun)%最小二乘拟合解函数例：已知一组测量数据以下，并给定一组拟合基函数y=1,y=x,y=x^2,y=cosx,y=e^x,y=sinx试求其最小二乘拟合函数p=0.38280.4070-0.3901-0.45980.07650.5653Xi00.511.522.53yi00.47940.84150.98150.91260.59850.1645第28页2025/4/1329实现流程图第29页function[a,b]=LZXEC(x,y)%离散试验数据点线性最小二乘拟合%试验数据点x坐标向量：X%试验数据点y坐标向量：Y%拟合一次项系数：a%拟合常数项：bif(length(x)==length(y))n=length(x);elsedisp('x和y维数不相等！');return;end%维数检验A=zeros(2,2);A(2,2)=n;B=zeros(2,1);fori=1:nA(1,1)=A(1,1)+x(i)*x(i);A(1,2)=A(1,2)+x(i);B(1,1)=B(1,1)+x(i)*y(i);B(2,1)=B(2,1)+y(i);endA(2,1)=A(1,2);s=A\B;a=s(1);b=s(2);2025/4/1330第30页例X12345y1.51.843.45.7>>x=1:5;>>y=[1.51.843.45.7];>>[a,b]=LZXEC(x,y)a=1.0000b=0.28002025/4/1331第31页例.回到引言实例，从散点图能够看出，纤维强度和拉伸倍数之间近似线性关系，故可选取线性函数为拟合函数建立正规方程组，其基函数为依据内积公式，可得正规方程组为2025/4/1332第32页解得残差平方和：拟合曲线与散点关系如右图:即为所求最小二乘解。故2025/4/1333第33页例：金属铜温度x与电阻y，线性拟合matlab程序2025/4/1334第34页2025/4/1335第35页

线性拟合在物理试验中应用十分广泛，比如弹性介质杨氏模量测量中应变与应力关系，电阻电路中电流与电压关系等。

有些物理量之间在一定范围内是线性关系，也可使用线性拟合方法，只是要注意其适用范围。

还有一个情况是量物理量之间并不存在线性关系，但经过适当变换后可转化为线性关系。2025/4/1336第36页

惯用线性变换

函数变换后函数2025/4/1337第37页2025/4/1338某化学反应，依据试验得到生成物浓度与时间关系以下表，求浓度与时间t最小二乘法时间t12345678浓度y*10^-34.006.408.008.809.229.509.709.86时间t910111213141516浓度y*10^-310.0010.2010.3210.4210.5010.5510.5810.60双曲线型：1/y=a+b/t指数型：y=a*e^(b/t)Lny=lna+b/t第38页2025/4/1339t=1:16;%时间y=[4.006.408.008.809.229.509.709.8610.00...10.2010.3210.4210.5010.5510.5810.60]*1e-3;%浓度plot(t,y,'*')%绘制散点图%1/y=a+b/twfun=ones(1,2);%权系数phifun=@(x)[ones(size(x));1./x];%拟合基函数[A1,b1,p1]=Least_square(wfun,phifun,t,1./y)%最小二乘拟合求解F1=t./(fliplr(p1)*phifun(1./t));%拟合值D1=y-F1;%误差%绝对值最大误差D1_max=max(abs(D1))%D1_max=norm(D1,inf)%均方误差D1_ME=norm(D1)%y=a*exp(b/t)==>lny=lna+b/t[A2,b2,p2]=Least_square(wfun,phifun,t,log(y))%最小二乘拟合求解F2=exp(p2*phifun(t));%拟合值D2=y-F2;%误差%绝对值最大误差D2_max=max(abs(D2))%D2_max=norm(D2,inf)%均方误差D2_ME=norm(D2)第39页2025/4/1340双曲线型：A1=16.00003.38073.38071.5843b1=1.0e+003*1.83290.5289p1=80.1745162.7225双曲线型误差：D1_max=5.6037e-004D1_ME=0.0012第40页2025/4/1341指数型：指数型误差：A2=16.00003.38073.38071.5843b2=-75.2639-16.8223p2=-4.4807-1.0567D2_max=2.7715e-004D2_ME=3.4101e-004第41页3.2多元线性拟合

影响一个物理量原因，很多情况下不止一个，为了得到描述物理规律近似函数，就必须采取多元线性拟合。

设物理量y随x1，x2，…，xk等k个物理量而改变，即：

为了寻求描述物理规律近似函数，经过试验测得n组数据(普通n>k)，利用拟合方法近似函数。2025/4/1342第42页测得n组试验数据以下：i123……nx1ix11x12x13……x1nx2ix21x22x23……x2nx3ix31x32x33……x3n.....xki.....xk1.....xk2.....xk3.....Xknyiy1y2y3……yn2025/4/1343第43页

设近似函数为

与一元线性拟合思绪相同，由n组试验数据代入上式，得到n个方程式组成k元线性方程组，用最小二乘原理确定函数系数A0,A1,…,Ak,使yi与Yi偏差平方和最小。2025/4/1344第44页偏差为即(i=1,2,…,n)令2025/4/1345第45页求极值2025/4/1346第46页得到2025/4/1347第47页化简整理后可得

将第一个方程式中A0提出，代入其它各式后，关于A0,A1,A2,…,Ak正规方程组则为2025/4/1348第48页2025/4/1349(r,s=1,2,…,k)第49页

物理学及各科学技术领域都普遍存在非线性函数关系，对此不能直接使用线性拟合方法，对这类问题能够采取不一样方法处理。方法一：变换为线性拟合有些非线性函数能够经过变量替换，转化成线性函数关系，然后对替换变量进行线性拟合，最终再还原为原始物理量。但不是全部函数关系都可经过变量替换而转化成线性函数关系。方法二：多项式拟合

任何一个连续函数，在一个比较小邻域内均可用多项式任意迫近。所以在物理学许多问题中，不论物理量直接存在何种函数关系，都可用多项式进行数据拟合。2025/4/1350第50页直线：多项式：有理函数：指数：3.3非线性拟合2025/4/1351第51页Log-linear:Gaussians:2025/4/1352第52页

设有n组试验数据xi，yi，（i=1,2,…,n），用k次多项式拟合，设拟合方程为：即：偏差为2025/4/1353第53页偏差平方和为：取极小值即2025/4/1354第54页整理得：即即得到了正则方程组2025/4/1355第55页MATLAB程序实现functionA=multifit(X,Y,m)%离散试验数据点多项式曲线拟合%试验数据点x坐标向量：X%试验数据点y坐标向量：Y%拟合多项式次数：m%拟合多项式系数向量：AN=length(X);M=length(Y);if(N~=M)disp('数据点坐标不匹配！');return;endc(1:(2*m+1))=0;b(1:(m+1))=0;。。。。。。forj=1:(2*m+1)%求出c和bfork=1:Nc(j)=c(j)+X(k)^(j-1);if(j<(m+2))b(j)=b(j)+Y(k)*X(k)^(j-1);endendendC(1,:)=c(1:(m+1));fors=2:(m+1)C(s,:)=c(s:(m+s));endA=C\b';%直接求解法求出拟合系数2025/4/1356第56页多项式拟合也可转化为多元拟合，只要令一元非线性拟合转化为多元线性拟合解正则方程组2025/4/1357第57页x1

2

3

4y4101826例用多项式拟合函数：解：设得即2025/4/1358第58页记系数矩阵为，则故正规方程组为解得2025/4/1359拟合曲线：第59页MATLAB程序解法：>>x=1:4;>>y=[4101826];>>A=multifit(x,y,2)A=-1.50004.90000.5000>>m=1:0.01:4;>>n=-1.5+4.9.*m+0.5.*m.^2;>>plot(x,y,'-*',m,n,'--r')>>P=polyfit(x,y,2)P=0.50004.9000-1.50002025/4/1360第60页4、matlab拟合函数线性拟合函数格式：p=linefit(x,y)说明：x，y输入同维数据向量p输出拟合多项式系数向量2025/4/1361第61页例：试验测量取得以下数据，求此物理规律近似表示式2025/4/1362第62页2025/4/1363第63页多项式拟合函数格式：p=polyfit(x,y,m)说明：x，y输入同维数据向量m拟合多项式次数p输出拟合多项式系数向量2025/4/1364第64页例：多项式拟合2025/4/1365第65页2025/4/1366第66页2025/4/1367第67页2025/4/1368第68页2025/4/1369第69页非线性函数拟合格式：a=lsqcurvefit(‘fun’,a0,x,y)说明：x，y输入同维数据向量fun拟合函数文件a输出拟合函数系数向量a0a初值2025/4/1370第70页例：已知利用以下数据，求a,b,c,d2025/4/1371第71页2025/4/1372第72页2025/4/1373第73页2025/4/1374第74页2025/4/1375第75页2025/4/1376多元非线性拟合例子：反应动力学中Hougen-Watson模型是非线性模型，其模型以下，y=(b1*x2-x3/b5)/(1+b2*x1+b3*x2+b4*x3,其中，y为反应速率，其三个决定原因为x1氢气，x2（n-戊烷）x3（异戊烷），下表是试验数据，试建其回归模型，求出未知参数b1,b2,b3,b4,b5氢气4702854