函数逼近与拟合法.pptx

单击此处编辑母版标题样式,*,*,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,第五讲 函数迫近与拟正当,第1页,内容提要,引言,函数迫近,傅里叶迫近,最小二乘法拟合,最小二乘法,多元线性拟合,非线性拟合,MATLAB,拟合函数,小结,2025/10/27 周一,2,第2页,例：考查某种纤维强度与其拉伸倍数关系,下表是实际测定,24,个纤维样品强度与对应拉伸倍数统计,:,纤维强度随拉伸倍数增加而增加。,1,、引言,2025/10/27 周一,3,第3页,24,个点大致分布在一条直线附近。,故可认为强度,y,与拉伸倍数,x,主要关系应为线性关系：,必须找到一个度量标准来衡量什么曲线最靠近全部数据点。,2025/10/27 周一,4,第4页,2025/10/27 周一,5,第5页,在一个包含有很多数据点区间内结构插值函数，必定使用高次多项式。而高次插值多项式是不稳定。,因为数据本身存在误差，利用插值方法得到插值多项式必定保留了全部测量误差，造成插值函数与物理规律差异较大。,试验数据拟合能够克服插值方法在处理这类问题中存在缺点。,对这么数据采取上一讲介绍插值方法近似求描述物理规律解析函数，必定存在以下缺点：,2025/10/27 周一,6,第6页,试验数据拟合基本思想：,使近似函数尽可能靠近数据点，而不要求近似函数一定经过全部数据点。,试验数据拟合能够在一定精度内找出反应物理量间客观函数关系解析式。假如试验数据存在误差，这种做法能够部分抵消原来数据中测量误差，从而使所得到拟合函数更加好地反应物理规律。,2025/10/27 周一,7,第7页,利用拟合能够处理两类物理问题：,物理规律已知，但描述物理规律解析式中一些系数未知，能够利用试验方法取得了物理量之间关系，经过拟合方法，求出这些系数近似值。,物理规律未知，利用试验方法取得了物理量之间关系，经过拟合方法，得到一个近似解析式，用于描述物理规律。,拟合函数尽可能靠近数据点怎样实现？,2025/10/27 周一,8,第8页,2,、,函数迫近,在区间,a,b,上已知一连续函数,f(x),，假如该函数表示式太过于复杂不利于进行计算机运算，就会利用一个简单函数去近似,f(x),，这就是函数迫近问题。,假如,f(x),表示式未知，只知道描述,f(x),一条曲线，这就是曲线拟合问题。,和插值问题不一样，迫近和拟合并不要求迫近函数在已知点上值一定等于原函数函数值，而是按照某种标准使得二者差值到达最小。,2025/10/27 周一,9,第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/10/27 周一,10,第10页,Chebyshev(,切比雪夫,),迫近,当一个连续函数定义在区间,-1,，,1,上时，能够展开成为切比雪夫级数。,2025/10/27 周一,11,实际应用中，依据精度要求来截取有限项数,第11页,function f=Chebyshev(y,k,x0),%,用切比雪夫多项式迫近已知函数,%,已知函数：,y,%,迫近已知函数所需项数：,k,%,迫近点,x,坐标：,x0,%,求得切比雪夫迫近多项式或在,x0,处,迫近值：,f,syms t;,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-t2),t,-1,1)/pi;,c(2)=2*int(subs(y,findsym(sym(y),sym(t)*T(2)/sqrt(1-t2),t,-1,1)/pi;,f=c(1)+c(2)*t;,for i=3:k+1,T(i)=2*t*T(i-1)-T(i-2);,c(i)=2*int(subs(y,findsym(sym(y),sym(t)*T(i)/sqrt(1-t2),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);,else,f=vpa(f,6);,end,end,end,2025/10/27 周一,12,第12页,2025/10/27 周一,13,例：用切比雪夫公式（取,6,项）迫近函数,1/（2-x），,并求当,x=0.5,时函数值,函数准确值为,1/（2-0.5,）,=0.6667,，可见迫近结果比较靠近,第13页,离散周期函数傅里叶迫近,function c=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);,end,c(m)=c(m)/N;,end,2025/10/27 周一,14,第14页,例,N,1,2,3,4,5,6,Y,0.8415,0.9093,0.1411,-0.7568,-0.9589,-0.2794,y=0.8415 0.9093 0.1411-0.7568-0.9589-0.2794;,c=DFF(y,6),c=,Columns 1 through 4,-0.0926-0.5003i -0.0260-0.0194i -0.0251+0.0000i -0.0260+0.0194i,Columns 5 through 6,-0.0926+0.5003i -0.0172-0.0000i,2025/10/27 周一,15,第15页,3.1,最小二乘法,首先，从一个简单例子来讨论一元线性拟合与最小二乘法问题。,为了含有普通性，把上式改写为：,经过试验测量，求金属铜电阻温度系数,，金属电阻与温度关系以下：,3,、最小二乘法拟合,2025/10/27 周一,16,第16页,经过试验测得金属铜温度,x,与电阻,y,数据以下：,x,i,(),Y,i,(,),x,i,(),Y,i,(,),x,i,(),Y,i,(,),0,4.38,70,5.58,140,6.74,10,4.56,80,5.74,150,6.94,20,4.70,90,5.96,160,7.12,30,4.86,100,6.06,170,7.28,40,5.08,110,6.26,180,7.42,50,5.24,120,6.44,190,7.60,60,5.40,130,6.58,200,7.78,2025/10/27 周一,17,第17页,设一元线性拟合函数为：,将试验数据代入,拟合函数，得到方程组,21,个线性方程,矛盾方程组,2025/10/27 周一,18,第18页,因为以上矛盾方程组不能确定一组唯一,A,0,和,A,1,，,也就是说，由方程组可求得,A,0,和,A,1,多组解，那么终究哪一组解最靠近客观真实值呢,?,按照拟合思想，应该使在每一个测量点拟合函数函数值尽可能靠近测量值，这么拟合函数才是满足要求，即：,定义偏差：,2025/10/27 周一,19,第19页,按照拟合思想，必须在每一个测量点偏差都很小，怎样到达这一要求？,不过因为偏差有正有负，求和时可能相互抵消，这并不能确保在每一个测量点偏差都很小。,方法一：偏差之和最小,尽管这种方法能够确保在每一个测量点偏差都很小，但这种方法数学处理比较困难。,方法二：偏差绝对值之和最小,2025/10/27 周一,20,第20页,这种方法既能够确保在每一个测量点偏差都很小，又方便数学处理，所以这种方法是可行。,方法三：偏差平方和最小,-,最小二乘法,2025/10/27 周一,21,第21页,残差向量各分量平方和记为：,最小二乘法：,以残差平方和最小问题解来确定拟合函数方法。,令,在回归分析中称为残差,(,i,=1,2,m,),残差向量：,2025/10/27 周一,22,第22页,由多元函数求极值必要条件，有,可得,即,2025/10/27 周一,23,第23页,上式为由,n,+1,个方程组成方程组，称正规方程组。,由,得,即,2025/10/27 周一,24,第24页,引入记号,则由内积概念可知,显然内积满足交换律,正规方程组便可化为,2025/10/27 周一,25,第25页,将其表示成矩阵形式：,其系数矩阵为对称阵。,所以正规方程组系数矩阵非奇异,即,依据,Crame,法则,正规方程组有唯一解，称其为最小二乘解。,2025/10/27 周一,26,第26页,2025/10/27 周一,27,function 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);,for i=1:n,for j=1:n,A(i,j)=sum(wfun(i).*phifun(i,:).*phifun(j,:);,end,b(i)=sum(wfun(i).*phifun(i,:).*y);,end,b=b;,a=Ab;p=a;,MATLAB,实现,第27页,2025/10/27 周一,28,x=0:0.5:3;%x,轴数据,y=0 0.4794 0.8415 0.9815 0.9126 0.5985 0.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)%,最小二乘拟合求解,syms x,digits(4)%,设定精度,Phifun=1;x;x.2;cos(x);exp(x);sin(x);,y=vpa(p*Phifun)%,最小二乘拟合解函数,例：已知一组测量数据以下，并给定一组拟合基函数,y=1,y=x,y=x2,y=cosx,y=ex,y=sinx,试求其最小二乘拟合函数,p=0.3828 0.4070 -0.3901 -0.4598 0.0765 0.5653,Xi,0,0.5,1,1.5,2,2.5,3,yi,0,0.4794,0.8415,0.9815,0.9126,0.5985,0.1645,第28页,2025/10/27 周一,29,实现流程图,第29页,function a,b=LZXEC(x,y),%,离散试验数据点线性,最小二乘拟合,%,试验数据点,x,坐标向量：,X,%,试验数据点,y,坐标向量：,Y,%,拟合一次项系数：,a,%,拟合常数项：,b,if(length(x)=length(y),n=length(x);,else,disp(x,和,y,维数不相等！,);,return;,end%,维数检验,A=zeros(2,2);,A(2,2)=n;,B=zeros(2,1);,for i=1:n,A(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);,end,A(2,1)=A(1,2);,s=AB;,a=s(1);,b=s(2);,2025/10/27 周一,30,第30页,例,X,1,2,3,4,5,y,1.5,1.8,4,3.4,5.7,x=1:5;,y=1.5 1.8 4 3.4 5.7;,a,b=LZXEC(x,y),a=1.0000,b=0.2800,2025/10/27 周一,31,第31页,例,.,回到引言实例,，从散点图能够看出，,纤维强度和拉伸倍数之间近似线性关系，故可选取线性函数,为拟合函数建立正规方程组，其基函数为,依据内积公式，可得,正规方程组为,2025/10/27 周一,32,第32页,解得,残差平方和：,拟合曲线与散点,关系如右图,:,即为所求最小二乘解。,故,2025/10/27 周一,33,第33页,例：金属铜温度,x,与电阻,y,，,线性拟合,matlab,程序,2025/10/27 周一,34,第34页,2025/10/27 周一,35,第35页,线性拟合在物理试验中应用十分广泛，比如弹性介质杨氏模量测量中应变与应力关系，电阻电路中电流与电压关系等。,有些物理量之间在一定范围内是线性关系，也可使用线性拟合方法，只是要注意其适用范围。,还有一个情况是量物理量之间并不存在线性关系，但经过适当变换后可转化为线性关系。,2025/10/27 周一,36,第36页,惯用线性变换,函数,变换后函数,2025/10/27 周一,37,第37页,2025/10/27 周一,38,某化学反应，依据试验得到生成物浓度与时间关系以下表，求浓度与时间,t,最小二乘法,时间,t,1,2,3,4,5,6,7,8,浓度,y*10-3,4.00,6.40,8.00,8.80,9.22,9.50,9.70,9.86,时间,t,9,10,11,12,13,14,15,16,浓度,y*10-3,10.00,10.20,10.32,10.42,10.50,10.55,10.58,10.60,双曲线型：,1/y=a+b/t,指数型：,y=a*e(b/t),Lny=lna+b/t,第38页,2025/10/27 周一,39,t=1:16;%,时间,y=4.00 6.40 8.00 8.80 9.22 9.50 9.70 9.86 10.00.,10.20 10.32 10.42 10.50 10.55 10.58 10.60*1e-3;%,浓度,plot(t,y,*)%,绘制散点图,%1/y=a+b/t,wfun=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/10/27 周一,40,双曲线型：,A1=,16.0000 3.3807,3.3807 1.5843,b1=,1.0e+003*,1.8329,0.5289,p1=,80.1745 162.7225,双曲线型误差：,D1_max=,5.6037e-004,D1_ME=,0.0012,第40页,2025/10/27 周一,41,指数型：,指数型误差：,A2=,16.0000 3.3807,3.3807 1.5843,b2=,-75.2639,-16.8223,p2=,-4.4807 -1.0567,D2_max=,2.7715e-004,D2_ME=,3.4101e-004,第41页,3.2,多元线性拟合,影响一个物理量原因，很多情况下不止一个，为了得到描述物理规律近似函数，就必须采取多元线性拟合。,设物理量,y,随,x,1,，,x,2,，,，,x,k,等,k,个物理量而改变，即：,为了寻求描述物理规律近似函数，经过试验测得,n,组数据,(,普通,nk),，利用拟合方法近似函数。,2025/10/27 周一,42,第42页,测得,n,组试验数据以下：,i,1,2,3,n,x,1i,x,11,x,12,x,13,x,1n,x,2i,x,21,x,22,x,23,x,2n,x,3i,x,31,x,32,x,33,x,3n,.,.,.,.,.,x,ki,.,.,.,.,.,x,k1,.,.,.,.,.,x,k2,.,.,.,.,.,x,k3,.,.,.,.,.,X,kn,y,i,y,1,y,2,y,3,y,n,2025/10/27 周一,43,第43页,设近似,函数,为,与一元线性拟合思绪相同，由n组试验数据代入,上,式，得到n个方程式组成k元线性方程组，用最小二乘原理确定函数系数A,0,A,1,A,k,使y,i,与Y,i,偏差平方和最小。,2025/10/27 周一,44,第44页,偏差为,即,(i=1,2,n),令,2025/10/27 周一,45,第45页,求极值,2025/10/27 周一,46,第46页,得到,2025/10/27 周一,47,第47页,化简整理后可得,将第一个方程式中A,0,提出，代入其它各式后，关于A,0,A,1,A,2,A,k,正规方程组则为,2025/10/27 周一,48,第48页,2025/10/27 周一,49,(r,s=1,2,k),第49页,物理学及各科学技术领域都普遍存在非线性函数关系，对此不能直接使用线性拟合方法，对这类问题能够采取不一样方法处理。,方法一：变换为线性拟合,有些非线性函数能够经过变量替换，转化成线性函数关系，然后对替换变量进行线性拟合，最终再还原为原始物理量。但不是全部函数关系都可经过变量替换而转化成线性函数关系。,方法二：多项式拟合,任何一个连续函数，在一个比较小邻域内均可用多项式任意迫近。所以在物理学许多问题中，不论物理量直接存在何种函数关系，都可用多项式进行数据拟合。,2025/10/27 周一,50,第50页,直线：,多项式：,有理函数：,指数：,3.3,非线性拟合,2025/10/27 周一,51,第51页,Log-linear:,Gaussians:,2025/10/27 周一,52,第52页,设有n组试验数据x,i,，y,i,，（i=1,2,n），用k次多项式拟合，设拟合方程为：,即：,偏差为,2025/10/27 周一,53,第53页,偏差平方和,为：,取极小值,即,2025/10/27 周一,54,第54页,整理得,：,即,即得到了正则方程组,2025/10/27 周一,55,第55页,MATLAB,程序实现,function A=multifit(X,Y,m),%,离散试验数据点多项式曲线拟合,%,试验数据点,x,坐标向量：,X,%,试验数据点,y,坐标向量：,Y,%,拟合多项式次数：,m,%,拟合多项式系数向量：,A,N=length(X);,M=length(Y);,if(N=M),disp(,数据点坐标不匹配！,);,return;,end,c(1:(2*m+1)=0;,b(1:(m+1)=0;,。,for j=1:(2*m+1)%,求出,c,和,b,for k=1:N,c(j)=c(j)+X(k)(j-1);,if(j x=1:4;,y=4 10 18 26;,A=multifit(x,y,2),A=-1.5000 4.9000 0.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.5000 4.9000 -1.5000,2025/10/27 周一,60,第60页,4,、,matlab,拟合函数,线性拟合函数,格式：,p=linefit(x,y),说明：,x,，,y,输入同维数据向量,p,输出拟合多项式系数向量,2025/10/27 周一,61,第61页,例：试验测量取得以下数据，求此物理规律近似表示式,2025/10/27 周一,62,第62页,2025/10/27 周一,63,第63页,多项式拟合函数,格式：,p=polyfit(x,y,m),说明：,x,，,y,输入同维数据向量,m,拟合多项式次数,p,输出拟合多项式系数 向量,2025/10/27 周一,64,第64页,例：多项式拟合,2025/10/27 周一,65,第65页,2025/10/27 周一,66,第66页,2025/10/27 周一,67,第67页,2025/10/27 周一,68,第68页,2025/10/27 周一,69,第69页,非线性函数拟合,格式：,a=lsqcurvefit(fun,a0,x,y),说明：,x,，,y,输入同维数据向量,fun,拟合函数文件,a,输出拟合函数系数 向量,a0 a,初值,2025/10/27 周一,70,第70页,例：已知,利用以下数据，求,a,b,c,d,2025/10/27 周一,71,第71页,2025/10/27 周一,72,第72页,2025/10/27 周一,73,第73页,2025/10/27 周一,74,第74页,2025/10/27 周一,75,第75页,2025/10/27 周一,76,多元非线性拟合例子：反应动力学中,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,285,470,470,470,100,100,470,100,100,100,285,285,n,戊烷,300,80,300,80,80,190,80,190,300,300,80,300,190,异,戊烷,10,10,120,120,10,10,65,65,54,120,120,10,120,速率,8.55,3.79,4.82,0.02,2.75,14.39,2.54,4.35,13,8.5,0.05,11.32,3.13,多元非线性拟合：,beta,r,t=nlinfit(X,Y,fun,beta0,options),beta:,返回拟合系数；,r,：残差；,X，Y,为输入数据；,beta0,为预设拟合系数初始值；,Options,为控制参数项,第76页,2025/10/27 周一,77,reactants=470 285 470 470 470 100 100 470 100 100 100 285 285;,300 80 300 80 80 190 80 190 300 300 80 300 190;,10 10 120 120 10 10 65 65 54 120 120 10 120;%,反应物,rate=8.55,3.79,4.82,0.02,2.75,14.39,2.54,4.35,13,8.5,0.05,11.32,3.13;%,反应速率,f=(b,x)(b(1)*x(:,2)-x(:,3)/b(5)./(1+b(2)*x(:,1)+b(3)*x(:,2)+b(4)*x(:,3);%,定义拟合函数,b0=1;0.1;0.2;0.1;2;%,初始值,betafit=nlinfit(reactants,rate,f,b0)%,回归系数求解,betafit=,1.2526,0.0628,0.0400,0.1124,1.1914,第77页,例：布氏硬度图像自动测量,2025/10/27 周一,78,第78页,布氏硬度测量关键是对压痕图像进行自动图像处理,以得到压痕直径,代入公式从而得到布氏硬度值。最好方法是抓取压痕边界,对边界用最小二乘法拟合出一个最贴近圆,求出该圆像素直径,从而求出真实直径。,2025/10/27 周一,79,第79页,2025/10/27 周一,80,第80页,2025/10/27 周一,81,第81页,2025/10/27 周一,82,第82页,小 结,函数迫近和拟合与插值法区分,线性最小二乘法应用和计算,最小二乘法程序设计,拟合函数了解与应用,2025/10/27 周一,83,第83页,谢 谢！,第84页,