4.插值法的程序.doc

1.1拉格朗日多项式和基函数的MATLAB程序function L=lagran1(X,Y)m=length(X); L=ones(m,m);for k=1: m V=1; for i=1:m if k=i V=conv(V,poly(X(i)/(X(k)-X(i); endendL1(k,:)=V; l(k,:)=poly2sym (V)endL=Y*l例1 给出节点数据， ，作五次拉格朗日插值多项式和基函数.解 1）首先将上述程序保存为M文件；2）然后在MATLAB工作窗口输入程序 X=-2.15 -1.00 0.01 1.02 2.03 3.25;Y=17.03 7.24 1.05 2.03 17.06 23.05;L= lagran1(X,Y)运行后输出五次拉格朗日插值多项式L及基函数ll = -0.0056*x5+0.0299*x4-0.0323*x3-0.0292*x2+0.0382*x-0.0004 0.0331*x5-0.1377*x4-0.0503*x3+0.6305*x2-0.4852*x+0.0048 -0.0693*x5+0.2184*x4+0.3961*x3-1.2116*x2-0.3166*x+1.0033 0.0687*x5-0.1469*x4-0.5398*x3+0.6528*x2+0.9673*x-0.0097 -0.0317*x5+0.0358*x4+0.2530*x3-0.0426*x2-0.2257*x+0.0023 0.0049*x5+0.0004 *x4-0.0266*x3+0.0001*x2+0.0220*x-0.0002L =1.0954-4.5745*x+3.3960*x2+2.1076*x3+0.0648*x4-0.2169*x51.2 拉格朗日插值及其误差估计的MATLAB程序function y,R=lagranzi(X,Y,x,M)n=length(X); m=length(x);for i=1:m z=x(i);s=0.0; for k=1:n p=1.0; q1=1.0; c1=1.0;for j=1:n if j=kp=p*(z-X(j)/(X(k)-X(j); end q1=abs(q1*(z-X(j);c1=c1*j; end s=p*Y(k)+s; end y(i)=s;endR=M*q1/c1;例2 已知，用拉格朗日插值及其误差估计的MATLAB主程序求的近似值，并估计其误差.解 1）首先将上述程序保存为M文件；2）然后在MATLAB工作窗口输入程序 x=2*pi/9; M=1; X=pi/6 ,pi/4, pi/3;Y=0.5,0.7071,0.8660; y,R=lagranzi(X,Y,x,M)运行后输出插值y及其误差限R为y = R =0.6434 8.8610e-004.2 牛顿插值法的MATLAB程序function y,R,A,C,L=newdscg(X,Y,x,M)n=length(X); m=length(x);for t=1:m z=x(t); A=zeros(n,n);A(:,1)=Y;s=0.0; p=1.0; q1=1.0; c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)- A(i-1,j-1)/(X(i)-X(i-j+1); end q1=abs(q1*(z-X(j-1);c1=c1*j; end C=A(n,n);q1=abs(q1*(z-X(n);for k=(n-1):-1:1C=conv(C,poly(X(k);d=length(C);C(d)=C(d)+A(k,k);end y(k)= polyval(C, z);endR=M*q1/c1;L(k,:)=poly2sym(C);例3 给出节点数据，作三阶牛顿插值多项式，计算，并估计其误差.解 1).首先将名为newdscg.m的程序保存为M文件;2).然后在MATLAB工作窗口输入程序 syms M,X=-4,0,1,2; Y =27,1,2,17; x=-2.345; y,R,P=newdscg(X,Y,x,M)运行后输出插值y及其误差限公式R，三阶牛顿插值多项式Py = 22.3211R =1323077530165133/562949953421312*M（即R =2.3503*M）P =11/12*x3+17/4*x2-25/6*x+13 分段线性插值的MATLAB程序function s= xxczhjt1(x0,y0,xi,x,y)s= interp1(x0,y0,xi); Sn= interp1(x0,y0,x0);plot(x0,y0,o,x0,Sn,-,xi,s,*,x,y,-.)legend(节点(xi,yi),分段线性插值函数Sn(x),插值点(x,s),被插值函数y)例4 设函数,在区间上取等距节点，构造分段线性插值函数,计算各小区间中点处的值，作出节点，插值点，和的图形；解 1).首先将名为newdscg.m的程序保存为M文件;2)然后在MATLAB工作窗口输入程序h=2*pi/9; x0=-pi:h:pi; y0=tan(cos(sqrt(3)+sin(2*x0)./(3+4*x0.2);xi=-pi+h/2:h:pi-h/2; fi=tan(cos(3(1/2)+sin(2*xi)./(3+4*xi.2);b=max(x0); a=min(x0); x=a:0.001:b; y=tan(cos(sqrt(3)+sin(2.*x)./(3+4*x.2);si=xxczhjt1(x0,y0,xi,x,y);Ri=abs(fi-si)./fi);i=xi,fi,si,Rititle(y=tan(cos(sqrt(3)+sin(2 x)/(3+4 x2)的分段线性插值的有关图形)运行后屏幕显示Ri （略），并且作出节点，插值点，和的图形(略)i = -2.7925 1.5492 1.5473 0.0012 -2.0944 1.5304 1.5330 0.0017 -1.3963 1.5294 1.5300 0.0004 -0.6981 1.5191 1.4687 0.0332 0 1.1110 1.1932 0.0740 0.6981 1.1455 1.1745 0.0253 1.3963 1.4962 1.4544 0.0279 2.0944 1.5544 1.5496 0.00312.7925 1.5557 1.5553 0.00024曲线拟合的线性最小二乘法及其MATLAB程序例6 给出一组数据点，试用线性最小二乘法求拟合曲线，并作出拟合曲线.xi-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6yi-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04解 （1）在MATLAB工作窗口输入程序 x=-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6; y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04;plot(x,y,r*),legend(实验数据(xi,yi)xlabel(x), ylabel(y),title(例7.2.1的数据点(xi,yi)的散点图)运行后屏幕显示数据的散点图（略）. （2）编写下列MATLAB程序计算在处的函数值，即输入程序 syms a1 a2 a3 a4x=-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6; fi=a1.*x.3+ a2.*x.2+ a3.*x+ a4运行后屏幕显示关于a1,a2, a3和a4的线性方程组fi = -125/8*a1+25/4*a2-5/2*a3+a4, -4913/1000*a1+289/100*a2-17/10*a3+a4, -1331/1000*a1+121/100*a2-11/10*a3+a4, -64/125*a1+16/25*a2-4/5*a3+a4, a4, 1/1000*a1+1/100*a2+1/10*a3+a4, 27/8*a1+9/4*a2+3/2*a3+a4, 19683/1000*a1+729/100*a2+27/10*a3+a4, 5832/125*a1+324/25*a2+18/5*a3+a4编写构造误差平方和的MATLAB程序 y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04;fi=-125/8*a1+25/4*a2-5/2*a3+a4, -4913/1000*a1+289/100*a2-17/10*a3+a4, -1331/1000*a1+121/100*a2-11/10*a3+a4, -64/125*a1+16/25*a2-4/5*a3+a4, a4, 1/1000*a1+1/100*a2+1/10*a3+a4, 27/8*a1+9/4*a2+3/2*a3+a4, 19683/1000*a1+729/100*a2+27/10*a3+a4, 5832/125*a1+324/25*a2+18/5*a3+a4;fy=fi-y; fy2=fy.2; J=sum(fy.2)运行后屏幕显示误差平方和如下J= (-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)2+(-4913/1000*a1+289/100*a2-17/10*a3+a4+171/2)2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)2+(a4+91/10)2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)2为求使达到最小，只需利用极值的必要条件 ，得到关于的线性方程组，这可以由下面的MATLAB程序完成，即输入程序 syms a1 a2 a3 a4J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)2+(-4913/1000*a1+289/100*a2-17/10*a3+a4.+171/2)2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)2+(a4+91/10)2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)2; Ja1=diff(J,a1); Ja2=diff(J,a2); Ja3=diff(J,a3); Ja4=diff(J,a4);Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3), Ja41=simple(Ja4),运行后屏幕显示J分别对a1, a2 ,a3 ,a4的偏导数如下Ja11=56918107/10000*a1+32097579/25000*a2+1377283/2500*a3+23667/250*a4-8442429/625Ja21 = 32097579/25000*a1+1377283/2500*a2+23667/250*a3+67*a4+767319/625Ja31 =1377283/2500*a1+23667/250*a2+67*a3+18/5*a4-232638/125Ja41 =23667/250*a1+67*a2+18/5*a3+18*a4+14859/25解线性方程组Ja11 =0，Ja21 =0，Ja31 =0，Ja41 =0，输入下列程序A=56918107/10000, 32097579/25000, 1377283/2500, 23667/250; 32097579/25000, 1377283/2500, 23667/250, 67; 1377283/2500, 23667/250, 67, 18/5; 23667/250, 67, 18/5, 18;B=8442429/625, -767319/625, 232638/125, -14859/25; C=B/A, f=poly2sym(C)运行后屏幕显示拟合函数f及其系数C如下C = 5.0911 -14.1905 6.4102 -8.2574f=716503695845759/140737488355328*x3-7988544102557579/562949953421312*x2+1804307491277693/2814749767106