MATLAB课件5-函数逼近与拟合法

科学计算与MATLAB主讲：唐建国中南大学材料科学与工程学院2013.09第五讲函数逼近与拟合法内容提要引言函数逼近傅里叶逼近最小二乘法拟合最小二乘法多元线性拟合非线性拟合MATLAB的拟合函数小结2020/1/113例：考察某种纤维的强度与其拉伸倍数的关系,下表是实际测定的24个纤维样品的强度与相应的拉伸倍数的记录:编号拉伸倍数强度编号拉伸倍数强度11.91.41355.5221.3145.2532.11.81565.542.52.5166.36.452.72.8176.5662.72.5187.15.373.531986.583.52.72087944218.98.51043.52298114.54.2239.58.1124.63.524108.1iiyxiiyx纤维强度随拉伸倍数增加而增加。1、引言2020/1/1141234567891012345678924个点大致分布在一条直线附近。故可认为强度y与拉伸倍数x的主要关系应为线性关系：xxy10)(为待定参数其中10,.),)(()(10越接近越好样本点与所有的数据点希望iiyxxxy必须找到一种度量标准来衡量什么曲线最接近所有数据点。2020/1/115012345678910110123456789101112131415数据点拟合曲线插值曲线YX2020/1/116①在一个包含有很多数据点的区间内构造插值函数，必然使用高次多项式。而高次插值多项式是不稳定的。②由于数据本身存在误差，利用插值方法得到的插值多项式必然保留了所有的测量误差，导致插值函数与物理规律差异较大。实验数据的拟合可以克服插值方法在处理这类问题中存在的缺点。对这样的数据采用上一讲介绍的插值方法近似求描述物理规律的解析函数，必然存在下列缺点：2020/1/117实验数据拟合的基本思想：使近似函数尽量靠近数据点，而不要求近似函数一定通过所有数据点。实验数据拟合可以在一定精度内找出反映物理量间客观函数关系的解析式。如果实验数据存在误差，这种做法可以部分抵消原来数据中的测量误差，从而使所得到的拟合函数更好地反映物理规律。2020/1/118利用拟合可以解决两类物理问题：1.物理规律已知，但描述物理规律的解析式中某些系数未知，可以利用实验方法获得了物理量之间的关系，通过拟合的方法，求出这些系数的近似值。2.物理规律未知，利用实验方法获得了物理量之间的关系，通过拟合的方法，得到一个近似的解析式，用于描述物理规律。拟合函数尽量靠近数据点如何实现？2020/1/1192、函数逼近在区间[a,b]上已知一连续函数f(x)，如果该函数表达式太过于复杂不利于进行计算机运算，就会利用一个简单函数去近似f(x)，这就是函数逼近问题。如果f(x)的表达式未知，只知道描述f(x)的一条曲线，这就是曲线拟合问题。和插值问题不同，逼近和拟合并不要求逼近函数在已知点上的值一定等于原函数的函数值，而是按照某种标准使得二者的差值达到最小。2020/1/1110逼近方法：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)2020/1/1111Chebyshev(切比雪夫)逼近当一个连续函数定义在区间[-1，1]上时，可以展开成为切比雪夫级数。00111()()()1,()()2()()nnnnnnfxfTxTxTxxTxxTxTx2020/1/1112functionf=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);endendend2020/1/1113离散周期函数的傅里叶逼近functionc=DFF(f,N)%用傅里叶级数逼近已知的离散周期函数%离散数据点：f%展开项数：N%离散傅里叶逼近系数：cc(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;end10nikxikyce2020/1/1114例N123456Y0.84150.90930.1411-0.7568-0.9589-0.2794y=[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.0000i2020/1/1115例：基于神经网络的高炉铁水硅含量预报模型根据RBF神经网络具有收敛速度快和全局优化的特点,建立了RBF网络模型,并将其应用对高炉铁水硅含量预报。监于铁水硅含量与炉缸温度之间的密切相关性,通过铁水硅含量来间接地反映炉内温度变化。采用MATLAB中的Newrbe函数进行函数逼近,对高炉一段连续时期内正常生产的数据的归一化处理后进行训练和仿真,提高了铁水硅含量预报的命中率。径向基函数RBF神经网络是一个只有一个隐藏层的三层前馈神经网络结构。它不仅具有很强的非线性映射能力,并且具有收敛速度快,全局优化的特点2020/1/11162020/1/11172020/1/11182020/1/11193.1最小二乘法首先，从一个简单的例子来讨论一元线性拟合与最小二乘法问题。0(1)yRx为了具有一般性，把上式改写为：0(1)RRt通过实验测量，求金属铜电阻温度系数α，金属电阻与温度关系如下：3、最小二乘法拟合2020/1/1120通过实验测得金属铜温度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.782020/1/1121设一元线性拟合函数为：0010101014.384.56104.70204.86307.78200AAAAAAAAA将实验数据代入拟合函数，得到方程组01YAAx0(1)yRx21个线性方程矛盾方程组2020/1/1122iy由于以上矛盾方程组不能确定一组唯一的A0和A1，也就是说，由方程组可求得A0和A1的多组解，那么究竟哪一组解最接近客观真实值呢?01iiYAAx按照拟合的思想，应当使在每一个测量点拟合函数的函数值尽量接近测量值，这样的拟合函数才是满足要求的，即：定义偏差：01iiiiiyYyAAx2020/1/11231()niiiyYx按照拟合的思想，必须在每一个测量点的偏差都很小，如何达到这一要求？但是由于偏差有正有负，求和时可能互相抵消，这并不能保证在每一个测量点的偏差都很小。方法一：偏差之和最小1()niiiyYx尽管这种方法可以保证在每一个测量点的偏差都很小，但这种方法数学处理比较困难。方法二：偏差绝对值之和最小2020/1/112421()niiiyYx这种方法既可以保证在每一个测量点的偏差都很小，又方便数学处理，所以这种方法是可行的。方法三：偏差的平方和最小-----最小二乘法2020/1/1125残差向量的各分量平方和记为：iiiyx)(minjiijjnyxaaaaS12010])([),,,(最小二乘法：以残差平方和最小问题的解来确定拟合函数的方法。令－－在回归分析中称为残差(i=1,2,…m)残差向量：222020/1/11260kaSminjiijjikyxax100])()[(2),,1,0(nk由多元函数求极值的必要条件，有可得即imiikmiijnjikjyxxxa110)()()(0])()()([10iikminjijikjyxxxa2020/1/1127上式为由n+1个方程组成的方程组，称正规方程组。njmiiikjmiijikyxaxx011)(])()([由得imiikmiijnjikjyxxxa110)()()(),,1,0(nk即miiikmiiniknmiiikmiiikyxxxaxxaxxa11111100)()()()()()()(),,1,0(nk2020/1/1128引入记号))(,),(),((21mrrrrxxx),,,(21myyyf),()(),(1ijmiikjkxx则由内积的概念可知imiikkyxf1)(),(),(),(kjjk显然内积满足交换律正规方程组便可化为),(),(),(),(1100faaaknknkk),,1,0(nk2020/1/1129将其表示成矩阵形式：的基，为函数类由于)(,),(),(10xxxn必然线性无关。因此)(,),(),(10xxxn其系数矩阵为对称阵。所以正规方程组的系数矩阵非奇异,即0])),det[((nnji根据Crame法则,正规方程组有唯一解，称其为最小二乘解。),(),(),(10fffn),(),(),(01000n),(),(),(11101n),(),(),(10nnnnnaaa102020/1/11302020/1/1131实现流程图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\