哈工大计算方法实验报告

ADDINCNKISM.UserStyle实验报告一题目（摘要）拉格朗日插值算法实验给定平面上个不同的数据点，，，；则满足条件，的次拉格朗日插值多项式是存在唯一的。若，且函数充足光滑，则当时，有误差预计式，前言：（目的和意义）实验目的：运用拉格朗日插值多项式求的近似值。输入：个数据点，；插值点输出：在插值点的近似值程序设计流程开始开始结束输入(xi,yi),ni=0,1,2,…,nL=0i=0xl=1xl=*xlj=0,1,…,i-1,i+1,…,nL=L+xl*i=n?输出y否是i=i+1实验成果、结论与讨论问题1拉格朗日插值多项式的次数越大越好吗？（1）设，N=5时，程序运行以下：TestLag(inline('1./(1+x.^2)'),-5,5,5,0.75:4.75);将区间[-5,5]分为了5段计算插值的点xi=0.75001.75002.75003.75004.7500计算出的插值yi=0.90540.52580.0096-0.3568-0.1595插值点处函数值yFact=0.64000.24620.11680.06640.0424计算误差err=-0.2654-0.27960.10720.42320.N=10时，程序运行以下：TestLag(inline('1./(1+x.^2)'),-5,5,10,0.75:4.75);将区间[-5,5]分为了10段计算插值的点xi=0.75001.75002.75003.75004.7500计算出的插值yi=0.69070.23300.11220.1084-0.2360插值点处函数值yFact=0.64000.24620.11680.06640.0424计算误差err=-0.05070.01320.0045-0.04200.2785N=20时，程序运行以下：TestLag(inline('1./(1+x.^2)'),-5,5,20,0.75:4.75);将区间[-5,5]分为了20段计算插值的点xi=0.75001.75002.75003.75004.7500计算出的插值yi=0.64130.24910.12820.19036.4150插值点处函数值yFact=0.64000.24620.11680.06640.0424计算误差err=-0.0013-0.0029-0.0114-0.1239-6.3726（2）设，N=5时，程序运行以下：TestLag(inline('exp(x)'),-1,1,5,[-0.95-0.050.050.95]);将区间[-1,1]分为了5段计算插值的点xi=-0.9500-0.05000.05000.9500计算出的插值yi=0.38630.95131.05122.5863插值点处函数值yFact=0.38670.95121.05132.5857计算误差err=1.0e-003*0.4471-0.10510.1069-0.6129N=10时，程序运行以下：TestLag(inline('exp(x)'),-1,1,10,[-0.95-0.050.050.95]);将区间[-1,1]分为了10段计算插值的点xi=-0.9500-0.05000.05000.9500计算出的插值yi=0.38670.95121.05132.5857插值点处函数值yFact=0.38670.95121.05132.5857计算误差err=1.0e-008*-0.3126-0.0055-0.0055-0.3714N=20时，程序运行以下：TestLag(inline('exp(x)'),-1,1,20,[-0.95-0.050.050.95]);将区间[-1,1]分为了20段计算插值的点xi=-0.9500-0.05000.05000.9500计算出的插值yi=0.38670.95121.05132.5857插值点处函数值yFact=0.38670.95121.05132.5857计算误差err=1.0e-012*0.73390-0.0002-0.5671问题2插值区间越小越好吗？（1）设，N=5时，程序运行以下：TestLag(inline('1./(1+x.^2)'),-1,1,5,[-0.95-0.050.050.95]);将区间[-1,1]分为了5段计算插值的点xi=-0.9500-0.05000.05000.9500计算出的插值yi=0.51360.99780.99780.5136插值点处函数值yFact=0.52560.99750.99750.5256计算误差err=0.0121-0.0002-0.00020.0121N=10时，程序运行以下：TestLag(inline('1./(1+x.^2)'),-1,1,10,[-0.95-0.050.050.95]);将区间[-1,1]分为了10段计算插值的点xi=-0.9500-0.05000.05000.9500计算出的插值yi=0.52430.99750.99750.5243插值点处函数值yFact=0.52560.99750.99750.5256计算误差err=0.00140.00000.00000.0014N=20时，程序运行以下：TestLag(inline('1./(1+x.^2)'),-1,1,20,[-0.95-0.050.050.95]);将区间[-1,1]分为了20段计算插值的点xi=-0.9500-0.05000.05000.9500计算出的插值yi=0.52560.99750.99750.5256插值点处函数值yFact=0.52560.99750.99750.5256计算误差err=1.0e-005*-0.70230.00000.0000-0.7023

（2）设，N=5时，程序运行以下：TestLag(inline('exp(x)'),-5,5,5,[-4.75-0.250.254.75]);将区间[-5,5]分为了5段计算插值的点xi=-4.7500-0.25000.25004.7500计算出的插值yi=-1.93211.42750.5882123.7146插值点处函数值yFact=0.00870.77881.2840115.5843计算误差err=1.9408-0.64870.6958-8.1303N=10时，程序运行以下：TestLag(inline('exp(x)'),-5,5,10,[-4.75-0.250.254.75]);将区间[-5,5]分为了10段计算插值的点xi=-4.7500-0.25000.25004.7500计算出的插值yi=0.04250.77961.2848115.6630插值点处函数值yFact=0.00870.77881.2840115.5843计算误差err=-0.0339-0.0008-0.0008-0.0788N=20时，程序运行以下：TestLag(inline('exp(x)'),-5,5,20,[-4.75-0.250.254.75]);将区间[-5,5]分为了20段计算插值的点xi=-4.7500-0.25000.25004.7500计算出的插值yi=0.00870.77881.2840115.5843插值点处函数值yFact=0.00870.77881.2840115.5843计算误差err=1.0e-007*-0.09140.00000.0000-0.1434问题3在区间考虑拉格朗日插值问题，为了使得插值误差较小，应如何选用插值节点？（1）设，N=5时，程序运行以下：TestLag2(inline('1./(1+x.^2)'),-1,1,5,[-0.95-0.050.050.95]);将区间[-1,1]分为了5段计算插值的点xi=-0.9500-0.05000.05000.9500计算出的插值yi=0.52540.99780.99780.5254插值点处函数值yFact=0.52560.99750.99750.5256计算误差err=1.0e-003*0.2071-0.3011-0.30110.2071N=10时，程序运行以下：TestLag2(inline('1./(1+x.^2)'),-1,1,10,[-0.95-0.050.050.95]);将区间[-1,1]分为了10段计算插值的点xi=-0.9500-0.05000.05000.9500计算出的插值yi=0.52550.99720.99720.5255插值点处函数值yFact=0.52560.99750.99750.5256计算误差err=1.0e-003*0.15620.26030.26030.1562N=20时，程序运行以下：TestLag2(inline('1./(1+x.^2)'),-1,1,20,[-0.95-0.050.050.95]);将区间[-1,1]分为了20段计算插值的点xi=-0.9500-0.05000.05000.9500计算出的插值yi=0.52560.99750.99750.5256插值点处函数值yFact=0.52560.99750.99750.5256计算误差err=1.0e-007*0.23180.23810.23810.2318（2）设，

N=5时，程序运行以下：TestLag2(inline('exp(x)'),-1,1,5,[-0.95-0.050.050.95]);将区间[-1,1]分为了5段计算插值的点xi=-0.9500-0.05000.05000.9500计算出的插值yi=0.38670.95141.05112.5857插值点处函数值yFact=0.38670.95121.05132.5857计算误差err=1.0e-003*0.0079-0.13170.1339-0.0108N=10时，程序运行以下：TestLag2(inline('exp(x)'),-1,1,10,[-0.95-0.050.050.95]);将区间[-1,1]分为了10段计算插值的点xi=-0.9500-0.05000.05000.9500计算出的插值yi=0.38670.95121.05132.5857插值点处函数值yFact=0.38670.95121.05132.5857计算误差err=1.0e-009*-0.5045-0.4791-0.4835-0.5994N=20时，程序运行以下：TestLag2(inline('exp(x)'),-1,1,20,[-0.95-0.050.050.95]);将区间[-1,1]分为了20段计算插值的点xi=-0.9500-0.05000.05000.9500计算出的插值yi=0.38670.95121.05132.5857插值点处函数值yFact=0.38670.95121.05132.5857计算误差err=1.0e-015*0.16650.3331-0.4441-0.8882

问题4考虑拉格朗日插值问题，内插比外推更可靠吗？程序运行以下：TestLag3([149],[550115185])计算插值的点xi=550115185计算出的插值yi=2.2667-20.2333-171.9000-492.7333插值点处函数值yFact=2.23617.071110.723813.6015计算误差err=-0.030627.3044182.6238506.3348程序运行以下：TestLag3([364964],[550115185])计算插值的点xi=550115185计算出的插值yi=3.11587.071810.167010.0388插值点处函数值yFact=2.23617.071110.723813.6015计算误差err=-0.8797-0.00070.55683.5626程序运行以下：TestLag3([100121144],[550115185])计算插值的点xi=550115185计算出的插值yi=4.43917.285010.722813.5357插值点处函数值yFact=2.23617.071110.723813.6015计算误差err=-2.2030-0.21390.00100.0658程序运行以下：TestLag3([169196225],[550115185])计算插值的点xi=550115185计算出的插值yi=5.49727.800110.800513.6006插值点处函数值yFact=2.23617.071110.723813.6015计算误差err=-3.2611-0.7291-0.07670.0009

实验所用函数functionyh=LagInterp(x,y,xh)%LagInterp计算拉格朗日插值%%Synopsis:yh=LagInterp(x,y,xh)%%Input:x=一维向量，将要做插值x的值%y=一维向量，将要做插值y的值%xh=数值或一维向量，计算插值的位置，支持计算一列xh的值%%Output:yh=数值或一维向量，通过计算插值的位置算出的插值ifmin(size(x))>1|min(size(y))>1 %判断x，y与否为一维向量error('x,ymustbevectors!');elseiflength(x)~=length(y) %判断x，y与否有同样多的元素error('xandymustagree!');endyh=zeros(size(xh));L=zeros(length(x)-1);forj=1:length(xh)fori=1:length(x)xCal=x;xCal(i)=[];%prod(xh(j)-xCal)/prod(x(i)-xCal)为拉格朗日基函数L(i)=prod(xh(j)-xCal)/prod(x(i)-xCal);yh(j)=yh(j)+L(i)*y(i); %yh=sum(L(i)*y(i))endend

functionTestLag(fx,a,b,n,xi)%TestLag实验题目11,2%%Synopsis:TestLag(fun,a,b,n,xi)%%Input:fx=用来验证插值计算精确率的函数%a,b=节点选用上下限%n=多项式次数，固定区间[-a,b]分段数%xi=要计算插值的点x=linspace(a,b,n);y=feval(fx,x);yi=LagInterp(x,y,xi);yFact=feval(fx,xi);err=yFact-yi;fprintf('将区间[%d,%d]分为了%d段\n',a,b,n);fprintf('计算插值的点xi=\n');disp(xi);fprintf('计算出的插值yi=\n');disp(yi);fprintf('插值点处函数值yFact=\n');disp(yFact);fprintf('计算误差err=\n');disp(err);

functionTestLag2(fx,a,b,n,xi)%TestLag2实验题目13%%Synopsis:TestLag2(fun,a,b,n,xi)%%Input:fx=用来验证插值计算精确率的函数%a,b=节点选用上下限%n=多项式次数，固定区间[-a,b]分段数%xi=要计算插值的点x=zeros(1,n);fork=1:nx(k)=cos((2*k-1)*pi/(2*n));%构造非等距节点endy=feval(fx,x);yi=LagInterp(x,y,xi);yFact=feval(fx,xi);err=yFact-yi;fprintf('将区间[%d,%d]分为了%d段\n',a,b,n);fprintf('计算插值的点xi=\n');disp(xi);fprintf('计算出的插值yi=\n');disp(yi);fprintf('插值点处函数值yFact=\n');disp(yFact);fprintf('计算误差err=\n');disp(err);

functionTestLag3(x,xi)%TestLag3实验题目14%%Synopsis:TestLag3(fun,a,b,n,xi)%%Input:x=构造Lagrange插值的节点%xi=要计算插值的点fx=inline('sqrt(x)');y=feval(fx,x);yi=LagInterp(x,y,xi);yFact=feval(fx,xi);err=yFact-yi;fprintf('计算插值的点xi=\n');disp(xi);fprintf('计算出的插值yi=\n');disp(yi);fprintf('插值点处函数值yFact=\n');disp(yFact);fprintf('计算误差err=\n');disp(err);

思考题：1.对实验1存在的问题，应如何解决？拉格朗日插值多项式的次数并不是越大越好，根据定义，插值式能够在节点处与实际函数匹配，但不能确保在节点之间较好的逼近实际函数。这个现象就是多项式摆动——Runge现象，有时多项式摆动能够通过谨慎选择基础函数的取样点来减小。普通采用分段插值能够较好的消除多项式摆动现象。2.对实验2存在的问题的回答，试加以阐明在分段段数相似的状况下，插值区间越大，误差越大。因素是大部分状况下，相对于比较大的区间，函数在比较小的区间上的函数值变化较缓和，因此即使出现摆动也不会偏离原函数太大。3.对实验3存在的问题的回答，试加以阐明第一问中已经提到多项式摆动能够通过谨慎选择基础函数的取样点来减小。我们来看下xk=cos(2k+1)π/2(n+1),k=0,1,…,n的分布和f(x)=1/(1+x2)的图像：能够看出，xk的分布是两端密集中间稀疏，f(x)的趋势是两边陡峭中间平缓，函数变化陡峭时节点增多正好能够增加插值的精确性。4.如何理解插值问题中的内插和外推？普通来说，内插时插值收敛于实际函数，一旦超出内插的范畴，插值函数会发散，且离插值区间越远外推误差越大。使用不用的插值办法在同一点外推的值也会相差诸多，这阐明外推本身就存在很大的不拟定性。实验报告二题目（摘要）龙贝格(Romberg)积分法考虑积分欲求其近似值，能够采用以下公式：（复化）梯形公式（复化）辛卜生公式（复化）柯特斯公式这里，梯形公式显得算法简朴，含有以下递推关系因此，很容易实现从低阶的计算成果推算出高阶的近似值，而只需要耗费较少的附加函数计算。但是，由于梯形公式收敛阶较低，收敛速度缓慢。因此，如何提高收敛速度，自然是人们极为关心的课题。为此，记为将区间进行等份的复化梯形积分成果，为将区间进行等份的复化辛卜生积分成果，为将区间进行等份的复化柯特斯积分成果。根据李查逊（Richardson）外推加速办法，可得到能够证明，如果充足光滑，则有（固定）这是一种收敛速度更快的一种数值求积公式，我们称为龙贝格积分法。该办法的计算可按下表进行……………很明显，龙贝格计算过程在计算机上实现时，只需开辟一种一维数组，即每次计算的成果，可寄存在位置上，其最后成果是寄存在位置上。前言：（目的和意义）运用龙贝格积分法计算积分程序设计流程1．准备初值，计算且（为等份次数）2．按梯形公式的递推关系，计算3．按龙贝格公式计算加速值4．精度控制。对给定的精度，若则终止计算，并取作为所求成果；否则，重复2~4步，直到满足精度为止。

实验成果、结论与讨论问题1运用龙贝格(Romberg)积分法计算积分程序运行以下：I=Romberginterg(inline('x.^2.*exp(x)'),0,1,25,1e-6)I=0.7183程序运行以下：I=Romberginterg(inline('exp(x).*sin(x)'),1,3,25,1e-6)I=10.9502程序运行以下：I=Romberginterg(inline('4./(1+x.^2)'),0,1,25,1e-6)I=3.1416程序运行以下：I=Romberginterg(inline('1./(1+x)'),0,1,25,1e-6)I=0.6931问题2被积函数无界，如何解决？程序运行以下：I=Romberginterg(inline('cos(x)./sqrt(1-x)'),0,1,25,1e-6)I=NaN由于函数在x=0处出现了0/0的状况，极限为1，因此Matlab的成果显示为NaN非数，解决办法是把下限0改为一种小数eps。I=Romberginterg(inline('sin(x)./x'),eps,1,25,1e-6)I=0.9461程序运行以下：I=Romberginterg(inline('cos(x)./sqrt(1-x)'),0,1,25,1e-6)I=NaN与（1）的因素相似，函数在x=1处出现了0/0的状况，成果为无穷，此时应选择变换t=x，将积分变为I=201I=20I=RombergInterg(inline('2*sin(x).*cos((sin(x)).^2)'),0,pi/2,25,1e-6)I=1.4996程序运行以下：I=Romberginterg(inline('cos(x)./sqrt(x)'),0,1,25,1e-6)I=NaN函数在x=0处出现了1/0的状况，成果为无穷。先将积分变为I=201cosxdxI=2*Romberginterg(inline('cos(x.^2)'),0,1,25,1e-6)I=1.8090程序运行以下：I=Romberginterg(inline('x.*sin(x)./(1-x.^2)'),-1,1,25,1e-6)I=NaN函数在x=1,-1处出现了sin1/0的状况，成果为无穷。被积函数为偶函数，做变换x=sint，积分变为I=2I=2*Romberginterg(inline('sin(x).*sin(sin(x))'),0,pi/2,25,1e-6)I=1.3825或者使用Gauss-Chebyshev求积公式：I=GaussInterg(inline('x.*sin(x)'),'Chebyshev',-1,1,1e-6)I=1.3825由于Gauss-Chebyshev求积公式求的是f(x)/1-x2在[-1,1]的积分，因此此处的f(x)问题3积分区间无限，如何解决？程序运行以下：I=Romberginterg(inline('exp(-x.^2)'),-10,10,25,1e-6)I=1.7725由于积分区间无限，而函数在[-10,10]外很小，能够用函数在[-10,10]上的积分近似这个无穷积分。或者使用Gauss-Hermite求积公式：I=GaussInterg(inline('x+1-x'),'Hermite',0,0,1e-6)I=1.7725由于Gauss-Hermite求积公式求的是e-x2f(x)在[-∞,∞]的积分，因此此处的程序运行以下：先做变换t=1/x，则原积分变为I=2*Romberginterg(inline('1./(1+x.^2)'),0,1,25,1e-6)I=1.5708程序运行以下：使用Gauss-Hermite求积公式：I=GaussInterg(inline('cos(x).^3'),'Hermite',0,0,1e-6)I=1.0820由于Gauss-Hermite求积公式求的是e-x2f(x)在[-∞,∞]的积分，因此此处的程序运行以下：使用Gauss-Laguerre求积公式：I=GaussInterg(inline('sin(x).^2'),'Laguerre',0,0,1e-6)I=0.4000由于Gauss-Laguerre求积公式求的是e-xf(x)在[0,∞]的积分，因此此处的f(x)为使用的函数functionI=RombergInterg(fun,a,b,npanel,tol,flag)%RombergInterg用Romberg办法求积分%%Synopsis:I=RombergInteg(fun,a,b,n,tol)%%Input:fun=(string)被积函数的函数名%a,b=积分下限和积分上限%npanel=(optional)将积分区间平分的段数，默认为25%tol=(optional)计算误差上限，默认为5e-9%flag=(optional)与否显示Romberg-T数表，默认为0为不显示%%Output:I=通过Romberg办法求积分的近似值ifnargin<4npanel=25;endifnargin<5tol=5e-9;endifnargin<6flag=0;endT(1,1)=TrapezoidInteg(fun,a,b,npanel);%T0(h)=T(h)err=1;%初始化误差值m=2;%初始化行计算的行序号whileerr>=tolT(m,1)=TrapezoidInteg(fun,a,b,2^m*npanel);%计算第m行T0(h/2^m)T(m,m)=0;%将矩阵T变成m*mforn=2:mT(m,n)=(4^(n-1)*T(m,n-1)-T(m-1,n-1))/(4^(n-1)-1);%Tm(h/2^k)与Tm-1(h/2^(k+1))和Tm-1(h/2^k)的递推关系enderr=abs(T(m,m)-T(m-1,m-1));%计算误差值m=m+1;%计算下一行endI=T(m-1,m-1);ifflag~=0disp(T);end

functionI=TrapezoidInteg(fun,a,b,npanel)%TrapezoidInteg用复化梯形公式求积分%%Synopsis:I=TrapezoidInteg(fun,a,b,n)%%Input:fun=(string)被积函数的函数名%a,b=积分下限和积分上限%npanel=(optional)将积分区间平分的段数，默认为25%%Output:I=通过复化梯形公式求积分的近似值ifnargin<4npanel=25;endnnode=npanel+1;%节点数=段数+1h=(b-a)/(nnode-1);%步长x=a:h:b;%将积分区间分段f=feval(fun,x); %求节点处被积函数的值I=h*(0.5*f(1)+sum(f(2:nnode-1))+0.5*f(nnode));

functionI=GaussInterg(fun,type,a,b,tol)%GaussInterg用Gauss型求积公式求积分,具体形式由使用者选用%%Synopsis:I=GaussInterg(fun,type,a,b)%I=GaussInterg(fun,type,a,b,tol)%%Input:fun=(string)被积函数的函数名%type=(string)具体Gauss求积公式的形式%a,b=积分上下限，Laguerre只计算0到inf，Hermite只计算-inf到inf，因此a,b对这两种形式无效%tol=(optional)误差容忍程度，默认为5e-5%%Output:I=通过Gauss型求积公式求积分的近似值ifnargin<5tol=5e-5;endn=7;%默认从7节点多项式开始计算IOld=1;%初始化IOld为1err=1;%初始化误差为1whileerr>=tolswitchtypecase'Legendre'%计算无奇点被积函数在-1到1的积分I=GaussLegendreInterg(fun,a,b,n);case'Chebyshev'%计算被积函数形如f(x)/sqrt(1-x^2)在-1到1的积分I=GaussChebyshevInterg(fun,a,b,n);case'Laguerre'%计算被积函数形如exp(-x)*f(x)的在0到inf的积分I=GaussLaguerreInterg(fun,n);case'Hermite'%计算被积函数形如exp(-x^2)*f(x)的在-inf到inf的积分I=GaussHermiteInterg(fun,n);otherwiseerror('Nosuchtype!');enderr=abs(I-IOld);%计算误差IOld=I;%把IOld赋值为I进行下次迭代n=n+1;%多项式节点递增end

functionI=GaussLegendreInterg(fun,a,b,n)%GaussLegendreInterg用Gauss-Legendre型求积公式计算无奇点被积函数在-1到1的积分%%Synopsis:I=GaussLegendreInterg(fun,a,b)%%Input:fun=(string)被积函数的函数名%a,b=积分下限和积分上限%n=(optional)高斯节点的个数，默认为7%%Output:I=通过Gauss型求积公式求积分的近似值ifnargin<4n=7;endifround(n)~=n|n<1error('Gauss节点个数必须为正整数');endp=LegendreIter(n);%构造n阶勒让德多项式p=p/polyval(p,1);%使Pn(1)=1;fori=1:n+1%构造p的导数p1p1(i)=(n-i+1)*p(i);endp1(n+1)=[];r=roots(p);%计算勒让德多项式的根L=polyval(p1,r);%计算与勒让德多项式的根匹配的高斯积分系数A=2./(1-r.^2)./L.^2;xi=((b-a)*r+(a+b))/2;%找出勒让德多项式的根在[a,b]上对应的数值xiyi=feval(fun,xi);%计算f(xi)I=(b-a)/2*(A'*yi);%I=sum(Ai*yi),前面的系数是积分区间变化是增加的常数

functionI=GaussChebyshevInterg(fun,a,b,n)%GaussChebyshevInterg用Gauss-Chebyshev型求积公式求积分计算被积函数形如f(x)/sqrt(1-x^2)在-1到1的积分%%Synopsis:I=GaussChebyshevInterg(fun,a,b)%%Input:fun=(string)被积函数的函数名%a,b=积分下限和积分上限%n=(optional)高斯节点的个数，默认为7%%Output:I=通过Gauss型求积公式求积分的近似值ifnargin<4n=7;endifround(n)~=n|n<1error('Gauss节点个数必须为正整数');endp=ChebyshevIter(n);%构造n阶切比雪夫多项式A=pi/n*ones(n,1);%计算与切比雪夫多项式的根匹配的高斯积分系数r=roots(p);%计算切比雪夫多项式的根xi=((b-a)*r+(a+b))/2;%找出切比雪夫多项式的根在[a,b]上对应的数值xiyi=feval(fun,xi);%计算f(xi)I=(b-a)/2*(A'*yi);%I=sum(Ai*yi),前面的系数是积分区间变化是增加的常数

functionI=GaussLaguerreInterg(fun,n)%GaussLaguerreInterg用Gauss-Laguerre型求积公式求积分计算被积函数形如exp(-x)*f(x)的在0到inf的积分%%Synopsis:I=GaussLaguerreInterg(fun)%I=GaussLaguerreInterg(fun,n)%%Input:fun=(string)被积函数的函数名%n=(optional)高斯节点的个数，默认为7%%Output:I=通过Gauss型求积公式求积分的近似值ifnargin<2n=7;endifround(n)~=n|n<1error('Gauss节点个数必须为正整数');endp=LaguerreIter(n);%构造n阶拉盖尔多项式p1=LaguerreIter(n+1);%构造n+1阶拉盖尔多项式xi=roots(p);%计算拉盖尔多项式的根L=polyval(p1,xi);%计算与拉盖尔多项式的根匹配的高斯积分系数A=xi./((n+1)^2*L.^2);yi=feval(fun,xi);%计算f(xi)I=A'*yi;%I=sum(Ai*yi),前面的系数是积分区间变化是增加的常数

functionI=GaussHermiteInterg(fun,n)%GaussHermiteInterg用Gauss-Hermite型求积公式求积分计算被积函数形如exp(-x^2)*f(x)的在-inf到inf的积分%%Synopsis:I=GaussHermiteInterg(fun)%I=GaussHermiteInterg(fun,n)%%Input:fun=(string)被积函数的函数名%n=(optional)高斯节点的个数，默认为7%%Output:I=通过Gauss型求积公式求积分的近似值ifnargin<2n=7;endifround(n)~=n|n<1error('Gauss节点个数必须为正整数');endp=HermiteIter(n);%构造n阶埃尔米特多项式p1=HermiteIter(n-1);%构造n-1阶埃尔米特多项式xi=roots(p);%计算埃尔米特多项式的根L=polyval(p1(1:n),xi);%计算与埃尔米特多项式的根匹配的高斯积分系数A=2^(n-1)*sqrt(pi)*factorial(n)./(n^2*L.^2);yi=feval(fun,xi);%计算f(xi)I=A'*yi;%I=sum(Ai*yi),前面的系数是积分区间变化是增加的常数

functionp=LegendreIter(n)%LegendreIter用递推的办法计算n次勒让德多项式的系数向量Pn+2(x)=(2*i+3)/(i+2)*x*Pn+1(x)-(i+1)/(i+2)*Pn(x)%%Synopsis:p=LegendreIter(n)%%Input:n=勒让德多项式的次数%%Output:p=n次勒让德多项式的系数向量ifround(n)~=n|n<0error('n必须是一种非负整数');endifn==0%P0(x)=1p=1;return;elseifn==1%P1(x)=xp=[10];return;endpBk=1;%初始化三项递推公式后项为P0pMid=[10];%初始化三项递推公式中项为P1fori=0:n-2pMidCal=zeros(1,i+3);%构造用于计算的x*Pn+1pMidCal(1:i+2)=pMid;pBkCal=zeros(1,i+3);%构造用于计算的PnpBkCal(3:i+3)=pBk;pFwd=(2*i+3)/(i+2)*pMidCal-(i+1)/(i+2)*pBkCal;%勒让德多项式三项递推公式Pn+2(x)=(2*i+3)/(i+2)*x*Pn+1(x)-(i+1)/(i+2)*Pn(x)pBk=pMid;%把中项变为后项进行下次迭代pMid=pFwd;%把前项变为中项进行下次迭代endp=pFwd;

functionp=ChebyshevIter(n)%ChebyshevIter用递推的办法计算n次勒让德-切比雪夫多项式的系数向量Tn+2(x)=2*x*Tn+1(x)-Tn(x)%%Synopsis:p=ChebyshevIter(n)%%Input:n=勒让德-切比雪夫多项式的次数%%Output:p=n次勒让德-切比雪夫多项式的系数向量ifround(n)~=n|n<0error('n必须是一种非负整数');endifn==0%T0(x)=1p=1;return;elseifn==1%T1(x)=xp=[10];return;endpBk=1;%初始化三项递推公式后项为T0pMid=[10];%初始化三项递推公式中项为T1fori=0:n-2pMidCal=zeros(1,i+3);%构造用于计算的x*Tn+1pMidCal(1:i+2)=pMid;pBkCal=zeros(1,i+3);%构造用于计算的PnpBkCal(3:i+3)=pBk;pFwd=2*pMidCal-pBkCal;%勒让德-切比雪夫多项式三项递推公式Tn+2(x)=2*x*Tn+1(x)-Tn(x)pBk=pMid;%把中项变为后项进行下次迭代pMid=pFwd;%把前项变为中项进行下次迭代endp=pFwd;

functionp=LaguerreIter(n)%LaguerreIter用递推的办法计算n次拉盖尔多项式的系数向量Ln+2(x)=(2*n+3-x)*Ln+1(x)-(n+1)*Ln(x)%%Synopsis:p=LaguerreIter(n)%%Input:n=拉盖尔多项式的次数%%Output:p=n次拉盖尔多项式的系数向量ifround(n)~=n|n<0error('n必须是一种非负整数');endifn==0%L0(x)=1p=1;return;elseifn==1%L1(x)=-x+1p=[-11];return;endpBk=1;%初始化三项递推公式后项为L0pMid=[-11];%初始化三项递推公式中项为L1fori=0:n-2pMidCal1=zeros(1,i+3);%构造用于计算的x*Ln+1(x)pMidCal1(1:i+2)=pMid;pMidCal2=zeros(1,i+3);%构造用于计算的Ln+1(x)pMidCal2(2:i+3)=pMid;pBkCal=zeros(1,i+3);%构造用于计算的Ln(x)pBkCal(3:i+3)=pBk;pFwd=((2*i+3)*pMidCal2-pMidCal1-(i+1)*pBkCal)/(i+2);%拉盖尔多项式三项递推公式Ln+2(x)=(2*n+3-x)*Ln+1(x)-(n+1)^2*Ln(x)pBk=pMid;%把中项变为后项进行下次迭代pMid=pFwd;%把前项变为中项进行下次迭代endp=pFwd;

functionp=HermiteIter(n)%HermiteIter用递推的办法计算n次埃尔米特多项式的系数向量Hn+2(x)=2*x*Hn+1(x)-2*(n+1)*Hn(x)%%Synopsis:p=HermiteIter(n)%%Input:n=埃尔米特多项式的次数%%Output:p=n次埃尔米特多项式的系数向量ifround(n)~=n|n<0error('n必须是一种非负整数');endifn==0%H0(x)=1p=1;return;elseifn==1%H1(x)=2*xp=[20];return;endpBk=1;%初始化三项递推公式后项为L0pMid=[20];%初始化三项递推公式中项为L1fori=0:n-2pMidCal=zeros(1,i+3);%构造用于计算的x*Hn+1(x)pMidCal(1:i+2)=pMid;pBkCal=zeros(1,i+3);%构造用于计算的Hn(x)pBkCal(3:i+3)=pBk;pFwd=2*pMidCal-2*(i+1)*pBkCal;%埃尔米特多项式三项递推公式Hn+2(x)=2*x*Hn+1(x)-2*(n+1)*Hn(x)pBk=pMid;%把中项变为后项进行下次迭代pMid=pFwd;%把前项变为中项进行下次迭代endp=pFwd;

思考题1.输入的参数有什么意义？输入的参数N越大，在有限区间上分段越小，由复化梯形公式的误差RTn=-((b-a))f''(η)/(12N2)，计算精度越高。2.在实验1中二分次数和精度的关系如何？二分次数越多，计算精度越高。由于每一次二分，N变为原来的二倍，由复化梯形公式的误差RTn=-((b-a))f''(η)/(12N2)，二分后误差约为二分前的四分之一。3.在实验2中给出的提示含有普遍性吗？存在其它的办法吗？试加以阐明。有普遍性，被积函数在积分区间内存在奇点。若在奇点上极限存在，则能够运用小数—机器Epsilon将奇点从积分区间去除，如问题2（1）。若奇点上极限不存在，则需要根据具体状况进行积分变换使被积函数在积分区间上不存在奇点。若求的是f(x)/1-x2在[-1,1]的积分，则能够4.在实验3中给出的提示含有普遍性吗？存在其它的办法吗？试加以阐明。有普遍性，积分区间均为无限区间。若求的是e-x2f(x)在[-∞,∞]的积分可用Gauss-Hermite求积公式，若求的是e-xf(x)在[0,∞]的积分可用Gauss-Laguerre实验报告三题目（摘要）牛顿(Newton)迭代法为初始猜想，则由递推关系产生逼近解的迭代序列，这个递推公式就是Newton法。当距较近时，很快收敛于。但当选择不当时，会造成发散。故我们事先规定迭代的最多次数。若超出这个次数，还不收敛，则停止迭代另选初值。前言：（目的和意义）运用牛顿迭代法求的根程序设计流程否否是否是定义输入开始输出迭代失败标志输出输出奇异标志结束

实验成果、结论与讨论问题1程序运行以下：r=NewtSolveOne('fun1_1',pi/4,1e-6,1e-4,10)r=0.7391程序运行以下：r=NewtSolveOne('fun1_2',0.6,1e-6,1e-4,10)r=0.5885问题2程序运行以下：r=NewtSolveOne('fun2_1',0.5,1e-6,1e-4,10)r=0.5671程序运行以下：r=NewtSolveOne('fun2_2',0.5,1e-6,1e-4,20)r=0.5669问题3程序运行以下：①p=LegendreIter(2)p=1.00000-0.3333p=LegendreIter(3)p=1.00000-0.60000p=LegendreIter(4)p=1.00000-0.857100.0857p=LegendreIter(5)p=1.00000-1.111100.23810② p=LegendreIter(6)p=1.00000-1.363600.45450-0.0216r=roots(p)'r=-0.150-0.2650.1530.264-0.23860.2386用二分法求根为：r=BinSolve('LegendreP6',-1,1,1e-6)r=-0.826-0.755-0.9590.0200.66160.959程序运行以下：①p=ChebyshevIter(2)p=1.00000-0.5000p=ChebyshevIter(3)p=1.00000-0.75000p=ChebyshevIter(4)p=1.00000-1.000000.1250p=ChebyshevIter(5)p=1.00000-1.250000.31250②p=ChebyshevIter(6)p=1.00000-1.500000.56250-0.0313r=roots(p)'r=-0.067-0.5480.0680.547-0.25880.2588用二分法求根为：r=BinSolve('ChebyshevT6',-1,1,1e-6)r=-0.163-0.755-0.8780.25880.0200.429与下列代码成果基本一致，只是元素次序稍有不同：j=0:5;x=cos((2*j+1)*pi/2/(5+1))x=0.0680.5480.2588-0.2588-0.547-0.068程序运行以下：①p=LaguerreIter(2)p=1-42p=LaguerreIter(3)p=1-918-6p=LaguerreIter(4)p=1-1672-9624p=LaguerreIter(5)p=1.0000-25.0000200.0000-600.0000600.0000-120.000②p=LaguerreIter(5)p=1.0000-25.0000200.0000-600.0000600.0000-120.000r=roots(p)'r=12.7327.8913.7111.45200.141用二分法求根为：r=BinSolve('LaguerreL5',0,13,1e-6)r=0.7221.47893.1507.72012.590程序运行以下：①p=HermiteIter(2)p=1.00000-0.5000p=HermiteIter(3)p=1.00000-1.50000p=HermiteIter(4)p=1.00000-3.000000.7500p=HermiteIter(5)p=1.00000-5.000003.75000②p=HermiteIter(6)p=1.00000-7.5000011.25000-1.8750r=roots(p)'r=-2.4872.488-1.6961.698-0.6170.616用二分法求根为：r=BinSolve('HermiteH6',-3,3,1e-6)r=-2.216-1.691-0.6030.8161.2442.104

所用到的函数functionr=NewtSolveOne(fun,x0,ftol,dftol,maxit)%NewtSolveOne用Newton法解方程f(x)=0在x0附近的一种根%%Synopsis:r=NewtSolveOne(fun,x0)%r=NewtSolveOne(fun,x0,ftol,dftol)%%Input:fun=(string)需规定根的函数及其导数%x0=猜想根，Newton法迭代初始值%ftol=(optional)误差，默认为5e-9%dftol=(optional)导数容忍最小值，不大于它表明Newton法失败，默认