函数逼近与拟合法讲义课件

第五讲函数逼近与拟合法内容提要引言函数逼近傅里叶逼近最小二乘法拟合最小二乘法多元线性拟合非线性拟合MATLAB的拟合函数小结2023/7/272例：考察某种纤维的强度与其拉伸倍数的关系,下表是实际测定的24个纤维样品的强度与相应的拉伸倍数的记录:纤维强度随拉伸倍数增加而增加。1、引言2023/7/27324个点大致分布在一条直线附近。故可认为强度y与拉伸倍数x的主要关系应为线性关系：必须找到一种度量标准来衡量什么曲线最接近所有数据点。2023/7/2742023/7/275在一个包含有很多数据点的区间内构造插值函数，必然使用高次多项式。而高次插值多项式是不稳定的。由于数据本身存在误差，利用插值方法得到的插值多项式必然保留了所有的测量误差，导致插值函数与物理规律差异较大。

实验数据的拟合可以克服插值方法在处理这类问题中存在的缺点。

对这样的数据采用上一讲介绍的插值方法近似求描述物理规律的解析函数，必然存在下列缺点：2023/7/276实验数据拟合的基本思想：

使近似函数尽量靠近数据点，而不要求近似函数一定通过所有数据点。

实验数据拟合可以在一定精度内找出反映物理量间客观函数关系的解析式。如果实验数据存在误差，这种做法可以部分抵消原来数据中的测量误差，从而使所得到的拟合函数更好地反映物理规律。2023/7/277利用拟合可以解决两类物理问题：物理规律已知，但描述物理规律的解析式中某些系数未知，可以利用实验方法获得了物理量之间的关系，通过拟合的方法，求出这些系数的近似值。物理规律未知，利用实验方法获得了物理量之间的关系，通过拟合的方法，得到一个近似的解析式，用于描述物理规律。拟合函数尽量靠近数据点如何实现？2023/7/2782、函数逼近

在区间[a,b]上已知一连续函数f(x)，如果该函数表达式太过于复杂不利于进行计算机运算，就会利用一个简单函数去近似f(x)，这就是函数逼近问题。如果f(x)的表达式未知，只知道描述f(x)的一条曲线，这就是曲线拟合问题。和插值问题不同，逼近和拟合并不要求逼近函数在已知点上的值一定等于原函数的函数值，而是按照某种标准使得二者的差值达到最小。2023/7/279逼近方法：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)2023/7/2710Chebyshev(切比雪夫)逼近当一个连续函数定义在区间[-1，1]上时，可以展开成为切比雪夫级数。2023/7/2711实际应用中，根据精度要求来截取有限项数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);endendend2023/7/27122023/7/2713例：用切比雪夫公式（取6项）逼近函数1/（2-x），并求当x=0.5时的函数值函数准确值为1/（2-0.5）=0.6667，可见逼近结果比较接近离散周期函数的傅里叶逼近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;end2023/7/2714例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.0000i2023/7/27153.1最小二乘法

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

通过实验测量，求金属铜电阻温度系数α，金属电阻与温度关系如下：3、最小二乘法拟合2023/7/2716

通过实验测得金属铜温度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.782023/7/2717

设一元线性拟合函数为：将实验数据代入拟合函数，得到方程组21个线性方程矛盾方程组2023/7/2718

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

按照拟合的思想，应当使在每一个测量点拟合函数的函数值尽量接近测量值，这样的拟合函数才是满足要求的，即：定义偏差：2023/7/2719

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

但是由于偏差有正有负，求和时可能互相抵消，这并不能保证在每一个测量点的偏差都很小。方法一：偏差之和最小

尽管这种方法可以保证在每一个测量点的偏差都很小，但这种方法数学处理比较困难。方法二：偏差绝对值之和最小2023/7/2720

这种方法既可以保证在每一个测量点的偏差都很小，又方便数学处理，所以这种方法是可行的。方法三：偏差的平方和最小-----最小二乘法2023/7/2721残差向量的各分量平方和记为：最小二乘法：以残差平方和最小问题的解来确定拟合函数的方法。令－－在回归分析中称为残差(i=1,2,…m)残差向量：2023/7/2722由多元函数求极值的必要条件，有可得即2023/7/2723上式为由n+1个方程组成的方程组，称正规方程组。由得即2023/7/2724引入记号则由内积的概念可知显然内积满足交换律正规方程组便可化为2023/7/2725将其表示成矩阵形式：其系数矩阵为对称阵。所以正规方程组的系数矩阵非奇异,即根据Crame法则,正规方程组有唯一解，称其为最小二乘解。2023/7/27262023/7/2727function[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实现2023/7/2728x=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.16452023/7/2729实现流程图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);2023/7/2730例X12345y1.51.843.45.7>>x=1:5;>>y=[1.51.843.45.7];>>[a,b]=LZXEC(x,y)a=1.0000b=0.28002023/7/2731例.回到引言的实例，从散点图可以看出，纤维强度和拉伸倍数之间近似线性关系，故可选取线性函数为拟合函数建立正规方程组，其基函数为根据内积公式，可得正规方程组为2023/7/2732解得残差平方和：拟合曲线与散点的关系如右图:即为所求的最小二乘解。故2023/7/2733例：金属铜温度x与电阻y，线性拟合matlab程序2023/7/27342023/7/2735

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

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

还有一种情况是量物理量之间并不存在线性关系，但经过适当变换后可转化为线性关系。2023/7/2736

常用的线性变换

函数变换后的函数2023/7/27372023/7/2738某化学反应，根据实验得到生成物的浓度与时间关系如下表，求浓度与时间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/t2023/7/2739t=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)2023/7/2740双曲线型：A1=16.00003.38073.38071.5843b1=1.0e+003*1.83290.5289p1=80.1745162.7225双曲线型误差：D1_max=5.6037e-004D1_ME=0.00122023/7/2741指数型：指数型误差：A2=16.00003.38073.38071.5843b2=-75.2639-16.8223p2=-4.4807-1.0567D2_max=2.7715e-004D2_ME=3.4101e-0043.2多元线性拟合

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

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

为了寻求描述物理规律的近似函数，通过实验测得n组数据(一般n>k)，利用拟合的方法近似函数。2023/7/2742测得n组实验数据如下：i123……nx1ix11x12x13……x1nx2ix21x22x23……x2nx3ix31x32x33……x3n.....xki.....xk1.....xk2.....xk3.....Xknyiy1y2y3……yn2023/7/2743

设近似函数为

与一元线性拟合的思路相同，由n组实验数据代入上式，得到n个方程式构成的k元线性方程组，用最小二乘原理确定函数系数A0,A1,…,Ak,使yi与Yi的偏差的平方和最小。2023/7/2744偏差为即(i=1,2,…,n)令2023/7/2745求极值2023/7/2746得到2023/7/2747化简整理后可得

将第一个方程式中的A0提出，代入其它各式后，关于A0,A1,A2,…,Ak的正规方程组则为2023/7/27482023/7/2749(r,s=1,2,…,k)

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

任何一个连续函数，在一个比较小的邻域内均可用多项式任意逼近。所以在物理学的许多问题中，不论物理量直接存在何种函数关系，都可用多项式进行数据拟合。2023/7/2750直线：多项式：有理函数：指数：3.3非线性拟合2023/7/2751Log-linear:Gaussians:2023/7/2752

设有n组实验数据xi，yi，（i=1,2,…,n），用k次多项式拟合，设拟合方程为：即：偏差为2023/7/2753偏差的平方和为：取极小值即2023/7/2754整理得：即即得到了正则方程组2023/7/2755MATLAB程序实现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';%直接求解法求出拟合系数2023/7/2756多项式拟合也可转化为多元拟合，只要令一元非线性拟合转化为多元线性拟合解正则方程组2023/7/2757x1

2

3

4y4101826例用多项式拟合函数：解：设得即2023/7/2758记系数矩阵为，则故正规方程组为解得2023/7/2759拟合曲线：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.50002023/7/27604、matlab拟合函数线性拟合函数格式：p=linefit(x,y)说明：x，y输入同维数据向量p输出拟合多项式的系数向量2023/7/2761例：实验测量获得如下数据，求此物理规律的近似表达式2023/7/27622023/7/2763多项式拟合函数格式：p=polyfit(x,y,m)说明：x，y输入同维数据向量m拟合多项式次数p输出拟合多项式的系数向量2023/7/2764例：多项式拟合2023/7/27652023/7/27662023/7/27672023/7/27682023/7/2769非线性函数拟合格式：a=lsqcurvefit(‘fun’,a0,x,y)说明：x，y输入同维数据向量fun拟合函数文件a输出拟合函数的系数向量a0a的初值2023/7/2770例：已知利用下列数据，求a,b,c,d2023/7/27712023/7/27722023/7/27732023/7/27742023/7/27752023/7/2776多元非线性拟合例子：反应动力学中Hougen-Watson的模型是非线性模型，其模型如下，y=(b1*x2-x3/b5)/(1+b2*x1+b3*x2+b4*x3,其中，y为反应速率，其三个决定因素为x1氢气，x2（n-戊烷）x3（异戊烷），下表是实验数据，试建其回归模型，求出未知参数b1,b2,b3,b4,b5氢气470