matlab实现插值法和曲线拟合

插值法和曲线拟合电子科技大学摘要：理解拉格朗日多项式插值、分段线性插值、牛顿前插，曲线拟合，用matlab编程求解函数，用插值法和分段线性插值求解同一函数，比较插值余项；用牛顿前插公式计算函数，计算函数值；对于曲线拟合，用不同曲线拟合数据。关键字：拉格朗日插值多项式；分段线性插值；牛顿前插；曲线拟合引言：在数学物理方程中，当给定数据是不同散点时，无法确定函数表达式，求解函数就需要很大的计算量，我们有多种方法对给定的表格函数进行求解，我们这里，利用插值法和曲线拟合对函数进行求解，进一步了解函数性质，两种方法各有利弊，适合我们进行不同的散点函数求解。正文：一、插值法和分段线性插值1拉格朗日多项式原理对某个多项式函数，已知有给定的k+1个取值点：其中对应着自变量的位置，而对应着函数在这个位置的取值。假设任意两个不同的xj都互不相同，那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为：其中每个为拉格朗日基本多项式（或称插值基函数），其表达式为：[3]拉格朗日基本多项式的特点是在上取值为1，在其它的点上取值为0。2分段线性插值原理给定区间[a,b],将其分割成a=x0x1…xn=b,已知函数y=f(x)在这些插值结点的函数值为yk=f(xk)(k=0,1,…,n)求一个分段函数Ih(x),使其满足：(1)Ih(xk)=yk，(k=0,1,…,n);(2)在每个区间[xk,xk+1]上,Ih(x)是个一次函数。易知,Ih(x)是个折线函数,在每个区间[xk,xk+1]上,(k=0,1,…,n)k1kk1k1kk1kkk,1)()()(xxxxxfxxxxxfxL，于是,Ih(x)在[a,b]上是连续的，但其一阶导数是不连续的。3拉格朗日插值多项式算法○1输入,(0,1,2,,)iixyin，令0)(xLn。○2对0,1,2,,in，计算0,()()/()nijijjjilxxxxx()()()nniiLxLxlxy4分段线性插值算法○1输入（xk,yk）,k=0,1,…,n;○2计算k1kk1k1kk1kkk,1)()()(xxxxxfxxxxxfxL5插值法和分段线性插值程序按下列数据分别作五次插值和分段线性插值，画出两条插值曲线以及给定数据点。求x1=0.32,x2=0.55,x3=0.68时的函数近似值，并比较两种方法的插值余项。ix0.300.420.500.580.660.72iy1.044031.084621.118031.156031.198171,23223拉格朗日插值程序：functionlagrintxi=[0.32,0.55,0.68];%xi=[0.2:0.001:0.8];x=[0.3,0.42,0.50,0.58,0.66,0.72];y=[1.04403,1.08462,1.11803,1.15603,1.19817,1.23223];L=zeros(size(y));m=length(xi);fori=1:mdxi=xi(i)-x;L(1)=prod(dxi(2:6))/prod(x(1)-x(2:6));L(6)=prod(dxi(1:6-1))/prod(x(6)-x(1:6-1));forj=2:6-1num=prod(dxi(1:j-1))*prod(dxi(j+1:6));den=prod(x(j)-x(1:j-1))*prod(x(j)-x(j+1:6));L(j)=num/den;endyi(i)=sum(y.*L);fprintf('x=%f,y=%f\n',xi(i),yi(i));endplot(xi,yi,'r');axis([0.20.81.031.24]);holdonplot(x,y,'b.','markersize',20)gridon分段线性插值算法程序：function[y]=div%xi=[0.3:0.001:0.72];x0=[0.3,0.42,0.50,0.58,0.66,0.72];y0=[1.04403,1.08462,1.11803,1.15603,1.19817,1.23223];k=1;xi=[0.32,0.55,0.68];forj=1:3fori=1:5ifxi(j)=x0(i)&&xi(j)=x0(i+1)&&k=3lx(1)=(xi(j)-x0(i+1))/(x0(i)-x0(i+1));lx(2)=(xi(j)-x0(i))/(x0(i+1)-x0(i));y(k)=lx(1)*y0(i)+lx(2)*y0(i+1);k=k+1;endendendplot(xi,y,'r');axis([0.20.81.031.24]);holdonplot(x0,y0,'b.','markersize',20)gridon6运算结果拉格朗日插值结果x=0.320000,y=1.049958x=0.550000,y=1.141271x=0.680000,y=1.209300拉格朗日插值余项：)72.0)(66.0)(58.0)(5.0)(42.0)(3.0(!6)()()()(655xxxxxxfxLxfxR分段插值结果ans=1.05081.14181.2095分段线性插值余项：)(max8))((2)(max)()(max11121,1xfhxxxxfxLxfRiiiiiixxxiiiixxxixxx由于拉格朗日插值的余项比分段线性插值的余项要求更为严格，点少、区间小的时候，拉格朗日插值要更好。但在区间较大、节点较多的时候，分段线性插值要更好。二、牛顿前插1牛顿前插原理n次牛顿前插公式：)1()1(!)1(!2!1)(!)(0020010000ntttnfttftffjtkfthxNnkjnkkn插值余项：)()!1()()1()()1(10nnnfhnntttthxR，),(0nxxm阶差分记作imimimfff111。[]ifxk阶差商是kmkmiiimmkxxxfxxxf0010)()(,,,差分和差商之间的关系是kikkiiihkfxxxf!,,,12牛顿前插算法○1输入),2,1,0(,,niyxnii。○2对nk,,3,2,1，ki,,2,1计算各阶差分imimimfff111○3计算函数值)1()1(!)1(!2!1)(!)(0020010000ntttnfttftffjtkfthxNnkjnkkn3牛顿前插程序：编写一个用牛顿前插公式计算函数值的程序，要求先输出差分表，再计算x点的函数值ix0.1250.2500.3750.5000.6250.750)(ixf0.7960.7730.7440.7040.6560.602分别求x=0.158和x=0.636的三次插值的值，并比较二者的插值余项。这里以x=0.636为例function[P]=newtonchax0=0.636;X=[0.1250.2500.3750.5000.6250.750];Y=[0.7960.7730.7440.7040.6560.602];h=abs(X(2)-X(1));n=find(abs(x0-X)3*h);X=X(n(1):n(end));Y=Y(n(1):n(end));w=length(X);R=zeros(w,w);R(:,1)=Y(:);fork=2:wforj=k:wR(j,k)=R(j,k-1)-R(j-1,k-1);endendt=(x0-X(1))/h;T=1;form=1:w-1T=T*(t-m+1);N(m)=R(m+1,m+1)*T/factorial(m);endP=R(1,1)+sum(N);4运行结果：差分表kxkfkfkf2kf3kf4kf50.796000000000000000000.773000000000000-0.02300000000000000000.744000000000000-0.029000000000000-0.0060000000000000000.704000000000000-0.040000000000000-0.011000000000000-0.005000000000000000.656000000000000-0.048000000000000-0.0080000000000000.0030000000000000.00800000000000000.602000000000000-0.054000000000000-0.0060000000000000.002000000000000-0.001000000000000-0.009000000000000X=0.636时ans=0.651459661824000x=0,158时ans=0.790229818880000三、曲线拟合1曲线拟合原理：给定数据(,),1,2,,jjxyjn。记拟合函数的形式为0011()()()()mmpxaxaxax（1.1），其中0()mkkx为已知的线性无关函数。求系数***01,,,maaa使得2201110(,,,)[()][()]nnmmjjkkjjjjkaaapxyaxy（1.2）取最小值。称**0()()mkkkpxax（1.3）为拟合函数或经验公式。如果),,1,0()(mkxxkk，则（1.3）为m次最小二乘拟合多项式2曲线拟合算法：已知数据对),,2,1)(,(njyxii，求多项式)()(0nmxaxPmiii，使得njjmiijimyxaaaa12010)(),,,(为最小。注意到此时kkxx)(，多项式系数maaa,,,10满足下面的线性方程组：mmmmmmmTTTaaaSSSSSSSSS10102112110其中)2,,2,1,0(1mkxSnjkjk，),,2,1,0(1mkxyTnjkjjk然后只要调用线性方程组的函数程序即可3曲线拟合程序：试分别用抛物线y=a+bx2和指数曲线y=aebx拟合下列数据ix12.53.54iy3.81.5026.033.0画出数据点和两条拟合曲线，并通过计算2个拟合函数残差向量的2范数来比较拟合优劣。用抛物线y=a+bx拟合程序：functionZXEx=[12.5^23.5^24^2];y=[3.81.5026.033.0];m=1;S=zeros(1,2*m+1);T=zeros(m+1,1);fork=1:2*m+1S(k)=sum(x.^(k-1));endfork=1:m+1T(k)=sum(x.^(k-1).*y);endA=zeros(m+1,m+1);a=zeros(m+1,1);fori=1:m+1forj=1:m+1A(i,j)=S(i+j-1);endenda=A\T;fork=1:m+1fprintf('a[%d]=%f\n',k,a(k));endp=polyfit(x,y,1);u=polyval(p,x);plot(sqrt(x),u,'b')holdonplot(sqrt(x),y,'b.')gridon指数曲线y=aebx拟合程序：functionZXE2x=[12.53.54];y=[3.81.5026.033.0];y=log(y);m=1;S=zeros(1,2*m+1);T=zeros(m+1,1);fork=1:2*m+1S(k)=sum(x.^(k-1));endfork=1:m+1T(k)=sum(x.^(k-1).*y);endA=zer