实验一函数插值方法报告

实验一函数插值方法报告一、问题提出对于给定的一元函数)(xfy的n+1个节点值(),0,1,,jjyfxjn。试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。数据如下：（1）jx0.40.550.650.800.951.05jy0.410750.578150.696750.901.001.25382求五次Lagrange多项式5L()x，和分段三次插值多项式，计算(0.596)f,(0.99)f的值。（提示：结果为(0.596)0.625732f,(0.99)1.05423f）（2）jx1234567jy0.3680.1350.0500.0180.0070.0020.001试构造Lagrange多项式6L()x，计算的(1.8)f,(6.15)f值。（提示：结果为(1.8)0.164762f,(6.15)0.001266f）二、要求1、利用Lagrange插值公式00,()nninkkiikkixxLxyxx编写出插值多项式程序；2、给出插值多项式或分段三次插值多项式的表达式；3、根据节点选取原则，对问题（2）用三点插值或二点插值，其结果如何；4、对此插值问题用Newton插值多项式其结果如何。Newton插值多项式如下：10010,()()[,,]()knnjkkjjkNxfxfxxxx其中：00,0()()[,,]kikiijjjikfxxxfxx三、目的和意义1、学会常用的插值方法，求函数的近似表达式，以解决其它实际问题；2、明确插值多项式和分段插值多项式各自的优缺点；3、熟悉插值方法的程序编制；4、如果绘出插值函数的曲线，观察其光滑性。四、实验学时：2学时五、实验步骤：1．进入C或matlab开发环境；2．根据实验内容和要求编写程序；3．调试程序；4．运行程序；5．撰写报告,讨论分析实验结果.解：一、编写插值函数结构程序Lagrange插值多项式M文件：lagrange1.mfunction[A1,LN,L1,B1]=lagrange1(X,Y)m=length(X);LN=ones(m,m);fork=1:mx1=1;fori=1:mifk~=ix1=conv(x1,poly(X(i)))/(X(k)-X(i));endendL1(k,:)=x1;B1(k,:)=poly2sym(x1)endA1=Y*L1;LN=Y*B1分段三次艾尔米特插值多项式的M文件：Hermite3.mfunction[f,ff]=Hermite3(x,y,y1)symst;f=0.0;if(length(x)==length(y))if(length(y)==length(y1))n=length(x);elsedisp('y和y的导数的维数不相等');return;endelsedisp('x和y的维数不相等！');return;endfori=1:nh=1.0;a=0.0;forj=1:nif(j~=i)h=h*(t-x(j))^2/((x(i)-x(j))^2);a=a+1/(x(i)-x(j));endendf=f+h*((x(i)-t)*(2*a*y(i)-y1(i))+y(i));endff=subs(f,'t');（1）、求五次Lagrange多项式5L()x，和分段三次插值多项式。在主显示区，输入五次Lagrange多项式5L()x程序：X=[0.40.550.650.800.951.05];Y=[0.410750.578150.696750.901.001.25382];[A1,LN,L1,B1]=lagrange1(X,Y)plot(X,A1);F=poly2sym(A1)运行后，输出五次Lagrange多项式5L()x的结果：A1=121.6264-422.7503572.5667-377.2549121.9718-15.0845F=(2139673480305281*x^5)/17592186044416-(1859275536318005*x^4)/4398046511104+(9836621836743*x^3)/17179869184-(414796119737013*x^2)/1099511627776+(2145751274873259*x)/17592186044416-1061478972867847/70368744177664拉格朗日插值多项式5L()x的图如下：在主显示区，输入分段三次艾尔米特插值多项式5L()x的程序：x=[0.40.550.650.800.951.05];y=[0.410750.578150.696750.901.001.25382];y1=[2.34400.90321.43290.99030.91705.1439];[f,ff]=Hermite3(x,y,y1);ff运行后，分段三次艾尔米特插值多项式5L()x的输出结果：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分段三次艾尔米特插值多项式5L()x的图如下：(2)、试构造Lagrange多项式6L()x。结果为(1.8)0.164762f,(6.15)0.001266f在主显示区，输入程序：X=[1234567];Y=[0.3680.1350.0500.0180.0070.0020.001];[A1,LN,L1,B1]=lagrange1(X,Y)plot(X,A1);F=poly2sym(A1)运行后，输出结果的Lagrange多项式6L()x的结果：A1=0.0001-0.00160.0186-0.11750.4419-0.96830.9950F=(4304240283865561*x^6)/73786976294838206464-(7417128346304051*x^5)/4611686018427387904+(223*x^4)/12000-(2821*x^3)/24000+(994976512675275*x^2)/2251799813685248-(19367*x)/20000+199/200Lagrange多项式6L()x的图如下：二、计算函数值计算函数值的主程序：lagrangezhi.mfunction[y,R]=lagrangezhi(X,Y,x,M)n=length(X);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;q1=1.0;c1=1.0;forj=1:nifj~=kp=p*(z-X(j))/(X(k)-X(j));endq1=abs(q1*(z-X(j)));c1=c1*j;ends=p*Y(k)+s;endy(i)=s;endR=M*q1/c1;（1）、计算(0.596)f、(0.99)f的值。在主显示区，输入程序：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）、计算(1.8)f、(6.15)f的值在主显示区，输入程序：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.0059在主显示区，输入程序：x=6.15;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.0013R=0.0042三、Newton插值多项式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;forj=2:nfori=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);fork=(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);symsMwcgs=M*Q/c1;Cw=q1/c1;在主显示区，输入5L()x的程序:x=[0.40.550.650.800.951.05];y=[0.410750.578150.696750.901.001.25382];[A,C,L,wcgs,Cw]=Newton(x,y)symsx;ezplot(L,[01.1]);运行结果如下，得到5L()xA=0.4108000000.57821.116000000.69671.18600.28000000.90001.35500.67600.9900001.00000.6667-2.2944-7.4261-15.302001.25382.53827.486124.451463.7551121.6264C=121.6264-422.7503572.5667-377.2549121.9718