数值分析实验报告

实验2.1多项式插值的振荡现象实验目的：在一个固定的区间上用插值逼近一个函数，显然Lagrange插值中使用的节点越多，插值多项式的次数就越高。我们自然关心插值多项式的次数增加时，Ln(x)是否也更加靠近被逼近的函数。Runge给出的一个例子是极著名并富有启发性的。实验内容：设区间[-1，1]上函数f(x)=1/(1+25x2)。考虑区间[-1,1]的一个等距划分，分点为xi=-1+2i/n，i=0，1，2，…，n，则拉格朗日插值多项式为201()()125nniiiLxlxx.其中，li(x)，i=0，1，2，…，n是n次Lagrange插值基函数。实验步骤与结果分析：实验源程序functionChap2Interpolation%数值实验二：“实验2.1：多项式插值的震荡现象”%输入：函数式选择，插值结点数%输出：拟合函数及原函数的图形promps={'请选择实验函数，若选f(x),请输入f,若选h(x),请输入h,若选g(x),请输入g:'};titles='charpt_2';result=inputdlg(promps,'charpt2',1,{'f'});Nb_f=char(result);if(Nb_f~='f'&Nb_f~='h'&Nb_f~='g')errordlg('实验函数选择错误！');return;endresult=inputdlg({'请输入插值结点数N:'},'charpt_2',1,{'10'});Nd=str2num(char(result));if(Nd1)errordlg('结点输入错误！');return;endswitchNb_fcase'f'f=inline('1./(1+25*x.^2)');a=-1;b=1;case'h'f=inline('x./(1+x.^4)');a=-5;b=5;case'g'f=inline('atan(x)');a=-5;b=5;endx0=linspace(a,b,Nd+1);y0=feval(f,x0);x=a:0.1:b;y=Lagrange(x0,y0,x);fplot(f,[ab],'co');holdon;plot(x,y,'b--');xlabel('x');ylabel('y=f(x)oandy=Ln(x)--');%--------------------------------------------------------------------functiony=Lagrange(x0,y0,x);n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nif(j~=k)p=p*(z-x0(j))/(x0(k)-x0(j));endends=s+p*y0(k);endy(i)=s;end实验结果分析(1)增大分点n=2，3，…时，拉格朗日插值函数曲线如图所示。n=6n=7n=8n=9n=10从图中可以看出，随着n的增大，拉格朗日插值函数在x=0附近较好地逼近了原来的函数f(x)，但是却在两端x=-1和x=1处出现了很大的振荡现象。并且，仔细分析图形，可以看出，当n为奇数时，虽然有振荡，但振荡的幅度不算太大，n为偶数时，其振荡幅度变得很大。通过思考分析，我认为，可能的原因是f(x)本身是偶函数，如果n为奇数，那么Lagrange插值函数Ln(x)的最高次项xn-1是偶次幂，比较符合f(x)本身是偶函数的性质；如果n为偶数，那么Lagrange插值函数Ln(x)的最高次项xn-1是奇次幂，与f(x)本身是偶函数的性质相反，因此振荡可能更剧烈。(2)将原来的f(x)换为其他函数如h(x)、g(x)，结果如图所示。其中h(x),g(x)均定义在[-5，5]区间上，h(x)=x/(1+x4)，g(x)=arctanx。h(x),n=7h(x),n=8h(x),n=9h(x),n=10g(x),n=7g(x),n=8g(x),n=9g(x),n=10分析两个函数的插值图形，可以看出：随着n的增大，拉格朗日插值函数在x=0附近较好地逼近了原来的函数f(x)，但是却在两端x=-5和x=5处出现了很大的振荡现象。并且，仔细分析图形，可以看出，当n为偶数时，虽然有振荡，但振荡的幅度不算太大，n为奇数时，其振荡幅度变得很大。原因和上面f(x)的插值类似，h(x)、g(x)本身是奇函数，如果n为偶数，那么Lagrange插值函数Ln(x)的最高次项xn-1是奇次幂，比较符合h(x)、g(x)本身是奇函数的性质；如果n为奇数，那么Lagrange插值函数Ln(x)的最高次项xn-1是偶次幂，与h(x)、g(x)本身是奇函数的性质相反，因此振荡可能更剧烈。实验3.1多项式最小二乘拟合实验目的：编制以函数{xk}k=0,…,n；为基的多项式最小二乘拟合程序。实验内容：对表中的数据作三次多项式最小二乘拟合。xi-1.0-0.50.00.51.01.52.0yi-4.447-0.4520.5510.048-0.4470.5494.552取权函数wi≡1，求拟合曲线**0nkkkx中的参数{αk}、平方误差δ2，并作离散据{xi,yi}的拟合函数的图形。实验源程序functionChap3CurveFitting%数值实验三:“实验3.1”%输出:原函数及求得的相应插值多项式的函数的图像以及参数alph和误差rx0=-1:0.5:2;y0=[-4.447-0.4520.5510.048-0.4470.5494.552];n=3;%n为拟合阶次alph=polyfit(x0,y0,n);y=polyval(alph,x0);r=(y0-y)*(y0-y)';%平方误差x=-1:0.01:2;y=polyval(alph,x);plot(x,y,'k--');xlabel('x');ylabel('y0*andpolyfit.y--');holdonplot(x0,y0,'*')gridon;disp(['平方误差:',num2str(r)])disp(['参数alph:',num2str(alph)])实验结果平方误差:2.1762e-005参数alph:1.9991-2.9977-3.9683e-0050.54912实验4.1实验目的：复化求积公式计算定积分.实验题目：数值计算下列各式右端定积分的近似值.实验要求：（1）若用复化梯形公式、复化Simpson公式和复化Gauss-LegendreI型公式做计算，要求绝对误差限为710*21，分别利用它们的余项对每种算法做出步长的事前估计.（2）分别用复化梯形公式，复化Simpson公式和复化Gauss-LegendreI型公式作计算.（3）将计算结果与精确解做比较，并比较各种算法的计算量.实验程序：1.事前估计的Matlab程序如下：（1）．用复化梯形公式进行事前估计的Matlab程序formatlonggx=2:0.01:3;f=-4*(3*x.^2+1)./(x.^2-1).^3;%二阶导函数%plot(x,f)%画出二阶导函数图像x=2.0;%计算导函数最大值f=-4*(3*x^2+1)/(x^2-1)^3;h2=0.5*10^(-7)*12/f;h=sqrt(abs(h2))%步长n=1/h;n=ceil(1/h)+1%选取的点数formatlonggx=0:0.01:1;f=8.*(3*x.^2-1)./(x.^2+1).^3;%二阶导函数%plot(x,f)%画出二阶导函数图像x=1;%计算导函数最大值f=8.*(3*x.^2-1)./(x.^2+1).^3;h2=0.5*10^(-7)*12/f;h=sqrt(abs(h2))%步长n=1/hn=ceil(1/h)+1%选取的点数formatlonggx=0:0.01:1;f=log(3).*log(3).*3.^x;%二阶导函数%plot(x,f);%画出二阶导函数图像x=1;%计算导函数最大值f=log(3)*log(3)*3^x;h2=0.5*10^(-7)*12/f;h=sqrt(abs(h2))%步长n=1/hn=ceil(1/h)+1%选取的点数formatlonggx=1:0.01:2;f=2.*exp(x)+x.*exp(x);%二阶导函数%plot(x,f)%画出二阶导函数图像x=2;%计算导函数最大值f=2.*exp(x)+x.*exp(x);h2=0.5*10^(-7)*12/f;h=sqrt(abs(h2))%步长n=1/hn=ceil(1/h)+1%选取的点数估计结果步长h及结点数n分别为h=0.000558156305651438n=1793h=0.000547722557505166n=1827h=0.000407071357304889n=2458h=0.000142479094906909n=7020（2）．用复化simpson公式进行事前估计的Matlab程序formatlonggx=2:0.01:3;f=-2*((-72*x.^2-24).*(x.^2-1)-192*x.^2.*(x.^2+1))./(x.^2-1).^5;%四阶导函数x=2.0;f=-2*((-72*x^2-24)*(x^2-1)-192*x^2*(x^2+1))/(x^2-1)^5;%计算导函数最大值h4=0.5*10^(-7)*180*16/f;h=sqrt(sqrt(abs(h4)))%步长n=1/h;%求分段区间个数n=2*ceil(1/h)+1%选取的点数formatlonggx=0:0.01:1;f=4*((-72*x.^2+24).*(x.^2+1)-192*x.^2.*(-x.^2+1))./(x.^2+1).^5;%四阶导函数x=1;f=4*((-72*x^2+24)*(x^2+1)-192*x^2*(-x^2+1))/(x^2+1)^5;%计算导函数最大值h4=0.5*10^(-7)*180*16/f;h=sqrt(sqrt(abs(h4)))%步长n=1/h;%求分段区间个数n=2*ceil(1/h)+1%选取的点数formatlonggx=0:0.01:1;f=log(3)^4*3.^x;%四阶导函数x=1;f=log(3)^4*3.^x;%计算导函数最大值h4=0.5*10^(-7)*180*16/f;h=sqrt(sqrt(abs(h4)))%步长n=1/h;%求分段区间个数n=2*ceil(1/h)+1%选取的点数formatlonggx=1:0.01:2;f=4*exp(x)+x.*exp(x);%四阶导函数plot(x,f)%画出原函数x=2;f=4*exp(x)+x.*exp(x);%计算导函数最大值h4=0.5*10^(-7)*180*16/f;h=sqrt(sqrt(abs(h4)))n=1/h;%求分段区间个数n=2*ceil(1/h)+1%选取的点数估计结果步长h及结点数n分别为h=0.0437490486013411n=47h=0.0588566191276542n=35h=0.0757645166218433n=29h=0.0424527247118546n=492.积分计算的Matlab程序：formatlonggpromps={'请选择积分公式，若用复化梯形，请输入T，用复化simpson，输入S，用复化Gauss_Legendre，输入GL：'};result=inputdlg(promps,'charpt4',1,{'T'});Nb=char(result);if(Nb~='T'&Nb~='S'&Nb~='GL')errordlg('积分公式选择错误');return;endresult=inputdlg({'请输入积分式题号1-4：'},'实验4.1',1,{'1'});Nb_f=str2num(char(result));if(Nb_f1|Nb_f4)errord