matlab实现插值法和曲线拟合

1、插值法和曲线拟合摘要：理解拉格朗日多项式插值、分段线性插值、牛顿前插，曲线拟合，用 matlab 编程求解函数，用插值法和分段线性插值求解同一函数，比较插值余项；用牛顿前插公式计算函数，计算函数值；对于曲线拟合，用不同曲线拟合数据。关键字：拉格朗日插值多项式；分段线性插值；牛顿前插；曲线拟合引言：正文：一、插值法和分段线性插值1 拉格朗日多项式原理对某个多项式函数，已知有给定的 k+ 1 个取值点：其中对应着自变量的位置，而对应着函数在这个位置的取值。其中每个其中对应着自变量的位置，而对应着函数在这个位置的取值。其中每个为拉格朗日基本多项式（或称插值基函数，其表达式为：3拉格朗日基本多项式的特

2、点是在拉格朗日基本多项式的特点是在上取值为1在其它的点0。给定区a,b, 将其分割成a=x0 x1=x0(i) & xi(j)=x0(i+1) & k=3 lx(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;end end endplot(xi,y,r);axis(0.2 0.8 1.03 1.24);hold on plot(x0,y0,b.,markersize,20) grid on运算结果拉格朗日插值结果x=0.320000

3、,y=1.049958 x=0.550000,y=1.141271 x=0.680000,y=1.209300拉格朗日插值余项：R5 (x) f (x) L5 (x)f 6 ( )(x 0.3)(x 0.42)(x 0.5)(x 0.58)(x 0.66)(x 0.72)6!分段插值结果ans =1.05081.14181.2095f (f ( i )2R max f ( x ) Lx x x1,i( x ) x x ( x xi )( x xi 1 )h2hi8 xi x xif ( x )由于拉格朗日插值的余项比分段线性插值的余项要求更为严格，点少、区间小的时候， 拉格朗日插值要更好。但在

4、区间较大、节点较多的时候，分段线性插值要更好。二、牛顿前插牛顿前插原理n 次牛顿前插公式：n0N n ( x0 th ) 0k k f k !k 1 ( t j ) f 0 j 0 f 0 t 1!2 2!0 t ( t 1) n n!0 t ( t 1) ( t n 1)插值余项：Rn ( x0 th ) t (t 1) (t n ) h n 1 f( n 1) () (x, xn )0( n1)!，0ii1m 阶差分记作m f m1 f ii1f xi k f x, x1 , , xk f ( xm )km 0 ( xm xi )i0imii差分和差商之间的关系是 ii, ,ik fk!h

5、k牛顿前插算法1n,xi , yi (i 0,2, n)。i2k ,2,3, ,n i ,2, ,k 计算各阶差分mfi m 1 fi m 1 fi3in0N n ( x0 th ) 0k k f k !k 1 ( t j ) f 0 j 0 f 0 t 1!2 2!0 t ( t 1) n n!0 t ( t 1) ( t n 1)牛顿前插程序：x点的函数值xi0.1250.2500.3750.5000.6250.750f ( xi )0.7960.7730.7440.7040.6560.602x=0.158x=0.636function P=newtoncha x0=0.636;X=0.1

6、25 0.250 0.375 0.500 0.625 0.750;Y=0.796 0.773 0.744 0.704 0.656 0.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);end endt=(x0-X(1)/h;T=1;for m=1:w-1 T=T*(t-m+1);N(m)=R(m+1,m+1)*T/factorial(m);

7、 endP=R(1,1)+sum(N);xkxkfkfk2 f k3 f k4 f k5fk0.796000000000000000000.77300000-0.023000000000000000000000000.74400000-0.02900000-0.006000000000000000000000000000000.70400000-0.04000000-0.01100000-0.005000000000000000000000000000000000000.65600000-0.04800000-0.008000000.0030000000.0080000000000000000

8、0000000000000000000000000.60200000-0.05400000-0.006000000.002000000-0.00100000-0.0090000000000000000000000000000000000000000000000X=0.636 时 ans =0.651459661824000 x=0,158 时 ans =0.790229818880000三、曲线拟合曲线拟合原理：给定数据( x j , y j ), j 1, 2, , n 。记拟合函数的形式为p(x) a00 (x) a1(x) amm (x)（1.1，k( x )其中k0 为已知的线性无关函

9、数。a * , a * , , a *求系数 01m 使得nnm ( a 0 , a1 , , am ) p ( x j ) y j 2 ak k ( x j ) y j 2j 1j 1 k 0（1.2）取最小值。称mp*(x) ak*k (xk 0（1.3）为拟合函数或经验公式。(x)x k(k如果 k，则（1.3）为次最小二乘拟合多项式曲线拟合算法：已知数据对( xi , y i )( j 1,2 , , n ) ，求多项式mP(x) ai xi(m ni 0，使得nm01m i jj ( a, a , ,) (a x i j 1 i 0) 2为最小。(x)x注意到此时 k，多项式系数 a

10、 0 , a1 , , am 满足下面的线性方程组： S 0S1S m a0 T0 SSS T 1 m 1 1 1 SSSS mm1 S 2 m aTm maT其中jjnnjjS k k jxj 1x(k 0,1,2,2m)， jy xk( k 0,1,2, , m)然后只要调用线性方程组的函数程序即可3 曲线拟合程序：xi12.53.54yi3.81.5026.033.0bxi12.53.54yi3.81.5026.033.0画出数据点和两条拟合曲线，并通过计算 2 个拟合函数残差向量的 2 范数来比较拟合优劣。用抛物线 y=a+bx 拟合程序： function ZXEx=1 2.52 3

11、.52 42;y=3.8 1.50 26.0 33.0; m=1;S=zeros(1,2*m+1);T=zeros(m+1,1);for k=1:2*m+1 S(k)=sum(x.(k-1); endfor k=1:m+1 T(k)=sum(x.(k-1).*y); endA=zeros(m+1,m+1);a=zeros(m+1,1); fori=1:m+1for j=1:m+1 A(i,j)=S(i+j-1);end end a=AT;for k=1:m+1 fprintf(a%d=%fn,k,a(k); endp=polyfit(x,y,1); u=polyval(p,x); plot(s

12、qrt(x),u,b)hold on plot(sqrt(x),y,b.) grid on指数曲线 y=aebx 拟合程序： function ZXE2x=1 2.5 3.5 4;y=3.8 1.50 26.0 33.0;y=log(y); m=1;S=zeros(1,2*m+1);T=zeros(m+1,1);for k=1:2*m+1 S(k)=sum(x.(k-1); endfor k=1:m+1 T(k)=sum(x.(k-1).*y); endA=zeros(m+1,m+1);a=zeros(m+1,1); fori=1:m+1for j=1:m+1 A(i,j)=S(i+j-1);end enda=AT;for k=1:m+1 fprintf(a%d=%fn,k,a(k); endp=polyfit(x,y,1); u=polyva