实验一函数插值方法.doc

数值分析课程设计实验报告实验一 函数插值方法一、问题提出对于给定的一元函数的n+1个节点值。试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。数据如下：0.40.550.650.800.951.050.410750.578150.696750.901.001.25382（1）求五次Lagrange多项式，和分段三次插值多项式，计算,的值。（提示：结果为,）12345670.3680.1350.0500.0180.0070.0020.001（2）试构造Lagrange多项式，计算的,值。（提示：结果为,）二、实验步骤1、利用Lagrange插值公式 编写出插值多项式程序； function A1,LN,L1,B1=lagrange(X,Y)m=length(X); LN=ones(m,m);for k=1: m x1=1; for i=1:m if k=i x1=conv(x1,poly(X(i)/(X(k)-X(i); endendL1(k,:)=x1; B1(k,:)=poly2sym (x1)endA1=Y*L1;LN=Y*B1在主显示区，输入五次Lagrange多项式L程序： X=0.4 0.55 0.65 0.80 0.95 1.05; Y=0.41075 0.57815 0.69675 0.90 1.00 1.25382; A1,LN,L1,B1=lagrange(X,Y) plot(X,A1); F=poly2sym(A1) 运行后，输出五次Lagrange多项式L的结果：A1 =121.6264 -422.7503 572.5667 -377.2549 121.9718 -15.0845F =(2139673480305281*x5)/17592186044416 - (1859275536318005*x4)/4398046511104 + (9836621836743*x3)- (414796119737013*x2)/1099511627776 + (2145751274873259*x)/17592186044416 - 1061478972867847/70368744177664拉格朗日插值多项式的图如下：2、给出插值多项式或分段三次插值多项式的表达式； function f,ff = Hermite3(x,y,y1)syms t;f = 0.0;if(length(x) = length(y) if(length(y) = length(y1) n = length(x); else disp(y和y的导数的维数不相等); return; endelse disp(x和y的维数不相等！ ); return;end for i=1:n h = 1.0; a = 0.0; for j=1:n if( j = i) h = h*(t-x(j)2/(x(i)-x(j)2); a = a + 1/(x(i)-x(j); end end f = f + h*(x(i)-t)*(2*a*y(i)-y1(i)+y(i);end ff = subs(f,t);在主显示区，输入分段三次艾尔米特插值多项式L的程序： x=0.4 0.55 0.65 0.80 0.95 1.05; y=0.41075 0.57815 0.69675 0.90 1.00 1.25382; y1=2.3440 0.9032 1.4329 0.9903 0.9170 5.1439; f,ff = Hermite3(x,y,y1); ff运行后，分段三次艾尔米特插值多项式L的输出结果：ff =(6400000000*(t - 4/5)2*(t - 11/20)2*(t - 13/20)2*(t - 19/20)2*(t - 21/20)2*(2240245151070481*t)/140737488355328 - 52393133567890089/8796093022208000)/184041 - (16000000*(6348013345609171*t)/140737488355328 - 85523418631741336287/1759218604441600000)*(t - 2/5)2*(t - 4/5)2*(t - 11/20)2*(t - 13/20)2*(t - 19/20)2)/169 + (16000000*(4105617466549689*t)/281474976710656 - 5238387122042657959/703687441776640000)*(t - 2/5)2*(t - 4/5)2*(t - 13/20)2*(t - 19/20)2*(t - 21/20)2)/9 - (256000000*(35097*t)/10000 - 46347/12500)*(t - 2/5)2*(t - 11/20)2*(t - 13/20)2*(t - 19/20)2*(t - 21/20)2)/81 - (400000000*(13147*t)/20000 - 449611/400000)*(t - 2/5)2*(t - 4/5)2*(t - 11/20)2*(t - 19/20)2*(t - 21/20)2)/81 - (10000000000*(84913*t)/11000 - 1833347/220000)*(t - 2/5)2*(t - 4/5)2*(t - 11/20)2*(t - 13/20)2*(t - 21/20)2)/9801分段三次艾尔米特插值多项式L的图如下：3、根据节点选取原则，对问题（2）用三点插值或二点插值，其结果如何； X=1 2 3 4 5 6 7; Y=0.368 0.135 0.050 0.018 0.007 0.002 0.001; A1,LN,L1,B1=lagrange(X,Y) plot(X,A1); F=poly2sym(A1)运行后，输出结果的Lagrange多项式L的结果：A1 =0.0001 -0.0016 0.0186 -0.1175 0.4419 -0.9683 0.9950F =(4304240283865561*x6)/73786976294838206464 - (7417128346304051*x5)/4611686018427387904 + (223*x4)/12000 - (2821*x3)/24000 + (994976512675275*x2)/2251799813685248 - (19367*x)/20000 + 199/200Lagrange多项式L的图如下：4、对此插值问题用Newton插值多项式其结果如何。Newton插值多项式如下： 其中： 计算函数值的主程序：lagrangezhi.mfunction y,R=lagrangezhi(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;（1）、计算 的值。在主显示区，输入程序： x=0.596; M=1; X=0.4,0.55,0.65,0.80,0.95,1.05; Y=0.41075,0.57815,0.69675,0.90,1.00,1.25382; y,R=lagrangezhi(X,Y,x,M)运行结果：y = 0.6257R = 2.2170e-008在主显示区，输入程序： x=0.99; M=1; X=0.4,0.55,0.65,0.80,0.95,1.05; Y=0.41075,0.57815,0.69675,0.90,1.00,1.25382; y,R=lagrangezhi(X,Y,x,M)运行结果：y = 1.0542R = 5.5901e-008（2）、计算f(1.8)的值 x=1.8; M=1; X=1,2,3,4,5,6,7; Y=0.368,0.135,0.050,0.018,0.007,0.002,0.001; y,R=lagrangezhi(X,Y,x,M)运行结果：y = 0.1648R = 0.0059Newton插值多项式Newton插值多项式主程序M文件：Newton.mfunction A,C,L,wcgs,Cw= Newton(X,Y)n=length(X); A=zeros(n,n); A(:,1)=Y;s=0.0; p=1.0; q=1.0; c1=1.0;for j=2:nfor i=j:nA(i,j)=(A(i,j-1)- A(i-1,j-1)/(X(i)-X(i-j+1);endb=poly(X(j-1);q1=conv(q,b); c1=c1*j; q=q1;endC=A(n,n); b=poly(X(n); q1=conv(q1,b);for k=(n-1):-1:1C=conv(C,poly(X(k); d=length(C); C(d)=C(d)+A(k,k);endL(k,:)=poly2sym(C); Q=poly2sym(q1);syms Mwcgs=M*Q/c1; Cw=q1/c1;在主显示区，输入的程序: x=0.4 0.55 0.65 0.80 0.95 1.05; y=0.41075 0.57815 0.69675 0.90 1.00 1.25382; A,C,L,wcgs,Cw= newploy(x,y) syms x; ezplot(L,0 1.1);运行结果如下，得到A = 0.4108 0 0 0 0 0 0.5782 1.1160 0 0 0 0 0.6967 1.1860 0.2800 0 0 0 0.9000 1.3550 0.6760 0.9900 0 0 1.0000 0.6667 -2.2944 -7.4261 -15.3020 0 1.2538 2.5382 7.4861 24.4514 63.7551 121.6264C = 121.6264 -422.7503 572.5667 -377.2549 121.9718 -15.0845L = (8558693921221117*x5)/70368744177664 - (3718551072636019*x4)/8796093022208 + (5036350380412441*x3)/8796093022208 - (3318368957896111*x2)/8796093022208 + (536437818718315*x)/4398046511104 - 8491831782942691/562949953421312 wcgs =(M*(x6 - (22*x5)/5 + (1583*x4)/200 - (3721*x3)/500 + (542206127247039*x2)/140737488355328 - (4682696525551953*x)/4503599627370496 + 4111390143022055/36028797018963968)/720Cw =0.0014 -0.0061 0.0110 -0.0103 0.0054 -0.0014 0.0002牛顿插值多项式的图如下：在主显示区，输入L的程序: x=1 2 3 4 5 6 7; y=0.368 0.135 0.050 0.018 0.007 0.002 0.001; A,C,L,wcgs,Cw= newploy(x,y) syms x; ezplot(L,0 8);运行结果如下，得到L：A = 0.3680 0 0 0 0 0 0 0.1350 -0.2330 0 0 0 0 0 0.0500 -0.0850 0.0740 0 0 0 0 0.0180 -0.0320 0.0265 -0.0158 0 0 0 0.0070 -0.0110 0.0105 -0.0053 0.0026 0 0 0.0020 -0.0050 0.0030 -0.0025 0.0007 -0.0004 0 0.0010 -0.0010 0.0020 -0.0003 0.0005 -0.0000 0.0001C =0.0001 -0.0016 0.0186 -0.1175 0.4419 -0.9683 0.9950L =(8608480567731121*x6)/147573952589676412928 - (7417128346304047*x5)/4611686018427387904 + (223*x4)/12000 - (2821*x3)/24000 + (7959812101402191*x2)/18014398509481984 - (19367*x)/20000 + 199/200wcgs =(M*(x7 - 28*x6 + 322*x5 - 1960*x4 + 6769*x3 - 13132*x2 + 13068*x - 5040)/5040Cw = 0.0002 -0.0056 0.0639 -0.38