2023年数值计算B大作业.pdf

珥流?充2普课 程 设 计课程名称：数值计算B设计题目：数值计算B大作业学 号：_姓 名：_完毕时间:题目一：多项式插值某气象观测站在8：0 0 (AM)开始每隔1 0 分钟对天气作如下观测，用三次多项式插值函数(N e wt on)逼近如下曲线，插值节点数据如上表，并求出9点 3 0 分该地区的温度(x=l 0 )。X12345678y2 2.52 3.32 4.42 1.702 5.22 8.524.82 5.4二、数学原理假设有n+1 个不同的节点及函数在节点上的值(x,y),(x,yn),插值多项式有如下形式：Pn(x)=aQ+a(x-x0)+a2(x-x0)(x-x,)+.+an(x-x0)(x-x,).(x-xn)(1 )其中系数/(i =0,1 ,2 n)为特定系数,可由插值样条P n(xp =y i (i=0,1,2 n)拟定。根据均差的定义，把 x 当作a,b上的一点，可得f(x)=f(x0)+f L xo，X (x-x0)fx,x0=f x0,xj+f x,X o，X,(X -X j )fx,x0,-xn.1 =f X ,X。，xn+f X ,x0,-xn(X X n)综合以上式子，把后一式代入前一式，可得到：f(x)=f x0+f x0，xj (x-xo)+f xo，X 1,x2(X -X()(X-X)+fx,X(p X (X-X 0)(X X I )+fx,X 0,X n ,x%+|(x)=N.(X )+4(x)其中N n (x)=f C xo+fxo，X (x-xo)+fx0，X,X,1 (X -XQ)(X-X j )+f E x,x0,-X(x-x0)(x-xn 4)(2)7?n(x)=f(X )N n(x)=fx,x0,X n ,X yn+l(x)(3)X (X -X,)+f x0,X,x2 (X -X,)(x-x2)+.+fE xo，X.X,(X -X j )(x-x2)-(x-x ).其中插值余项为：Rnn(x)=f(x)-N (x)=(n+1)!n+l(x)J介于X o，X,.X k之间。三、程序设计fu n c t i on y ,A,C,L=ne wd s c g(X,Y,x,M)%y为相应x的值,A为差商表，C为多项式系数，L为多项式%X为给定节点，Y为节点值，x为待求节点n=l e n g t h(X);m=le n g t h(x)；/n 为 X 的长度for t=l:mz=x(t)；A =z e r o s(n,n);A (:,1)=Y,;s=0 .0；p =1.0;ql=l.0;cl=l.0；f or j =2 :nfor i =j：nA(i,j )=(A(i,j -1)-A(i-1,j-l)/(X(i )-X(i -j +1);en dql=abs(q 1*(z -X(j-1 )；c 1=c l*j;e n dC=A(n,n);q l=ab s(ql*(z -X(n);f o r k=(n 1)：-1:1C=conv (C,p o 1 y (X (k);d=l e n g t h(C);C(d)=C(d)+A(k,k);y (k)=p o ly v al(C,z);%输出 y 值L (k,:)=p oly 2 s y m(C);为输出多项式 s y m s M,X=1,3,5,7 1 ;Y=2 2.5,2 4.4,2 5 .2,2 4.8；x=1 0;y,A,C,L=n e wd s c g(X,Y,x,M)2 1.7 3 1 32 2.5 0 0 02 4.4 0 0 00.9 5 0 02 5.2 0 2 30 .4 0 0 0 -0 .1 3 7 52 4.8 0 0 02 0 2 3-0.1 5 0 00 0 2 1c-0.00 2 1 -0.1 1 87 1.4 52 1 21.168 8L=-x-3/480-(1 9*x-2)/1 6 0+(69 7*x)/4 80+3387/1 6 0四、结果分析和讨论对于不超过三次的插值多项式,x假如选取1,3,5,7 这三个点可以得到较好的三次插值多项式 L=-0.0 0 2 lx-3-0.1 1 87x-2+l.4521x+2 1.1688。当 x=10 时，也即 9 点 30 分时的温度为21.7 317度，结果分析知此值应是偏小的。对于选取不同的插值节点，可以得到不同的插值多项式,误差也不尽相同。五、完毕题目的体会与收获牛顿插值法的重要一点就是对插值节点的选取，通过本题的编程很好的加深了对其概念的理解以及提高了应用牛顿插值法的能力，学会了运用Matl a b 软件对牛顿插值法相关问题进行编程求解，对 Matla b 计算方法与程序编辑更加熟悉。使我对这类问题的理解有了很大的提高。题目二：曲线拟合在某钢铁厂冶炼过程中，根据记录数据的含碳量与时间关系,试用最小二乘法拟合含碳量与时间力的拟合曲线，并绘制曲线拟合图形。t(分)0 5 10 1 5 20 25 3 0 3 5 40 4550 55(xKT4)0 1.27 2.16 2.86 3.4 4 3.87 4.1 5 4.37 4.51 4.58 4.02 4.64二 数学原理从整体上考虑近似函数P（x）同所给数据点（巧，）（i=0,1,m）误差，=PU,）-Z （i=0,1,m）的大小，常用的方法有以下三种:一是误差=以七）一（i=0,1,，m）绝对值mX 11的最大值盥N d,即误差向量，=，小*的8 范数；二是误差绝对值的和金“，即误m差向量r 的 1一范数;三是误差平方和，=。的算术平方根，即误差向量r 的2一范数;前两种方法简朴、自然，但不便于微分运算，后一种方法相称于考虑2 范数的平方，因此在曲线拟合mYr/中常采用误差平方和M 来度量误差。（i=0,1,，m）的整体大小。数据拟合的具体作法是:对给定数据（X-）（i=0,1,，0 1）,在取定的函数类中，求2 3）6%使误差”。（七）一必仕=0,1,1）的平方和最小，即由 2 （巧）-乂=一/=0=/=0从几何意义上讲,就是寻求与给定点（X，-）（i=0,1,，m）的距离平方和为最小的曲线y=（x）。函数P（x）称为拟合函数或最小二乘解，求拟合函数P（x）的方法称为曲线拟合的最小二乘法。在曲线拟合中，函数类中可有不同的选取方法。三、程序设计 x=0,5,1 0,15,20,25,3 0,35,40,4 5,5 0,55；y=0,1.2 7,2.16,2.8 6,3.44,3.8 7,4.1 5,4.3 7,4.51,4.5 8,4.02,4.6 4 *1 0 （-4）;p=p o ly fit(x,y,2)P1.Oe-0 4*-0.00240.20370.2305 p l o t (x,y,r)四、结果分析和讨论通过最小二乘法得到的曲线拟合多项式是P=(-0.0 024x-2+0.2 0 3 7x+0.2 3 05)*10、运用M a t 1 a b软件调用最小二乘法函数即可得到多项式拟合函数，由于数据较少得到的拟合曲线不够光滑,也许会存在一定的误差。拟合曲线总体呈现增函数趋势，与数据较为吻合。五、完毕题目的体会与收获曲线拟合较常用的方法之一就涉及最小二乘法，因此可以纯熟使用Ma tla b 进行数据的曲线拟合变得尤为重要。通过完毕本题的拟合过程,对于最小二乘法曲线拟合的操作更加的纯熟，可以较好的完毕各类数据的拟合。题目三:线性方程组求解分别运用LU分解;平方根法与改善平方根法;追赶法求解下列几个不同类型的线性方程组，并与准确值比较结果，分析误差产生因素。（1）设线性方程组42-3-1210000玉586-5-36501001242-2-132-1031不30-215-13-11942-426-167-3323X5386-8571726-35X64602-13-425301X7131610-11-917342-122/38462-7139201241900-18-3-24 863-1产。一-21/=（1,-1,0,1,2,0,3,1,-1,2/二、数学原理将A分解为一个下三角矩阵L和一个上三角矩阵U,即：A =L U,1 0 0一U U2,一 Wl/i其中 I 尸h 1 -0,u=“22，,W2M：1 l.00 o ;Jnl 4,2 L_ 0 0 ”解A x=b的问题就等价于规定解两个三角形方程组:右 Ly=b,求y；A Ux=y,求Xo设A为非奇异矩阵,且有分解式A=LU,L为单位下三角阵,U为上三角阵。L,U的元素可以有n步直接计算定出。用直接三角分解法解Ax=b （规定A的所有顺序主子式都不为零）的计算公式：ut i=,.（/=1,2,）,=%JU,i=2,3,，n.A 计算 U 的第 r 行，L 的第 r 列r-元素（i=2,3,，n）：urj=ari-lrkuki,i =r,r+1,n；k=/力=（%一/欢 以r）/rr，i=+1 ,，n,且 J T#!.k=求解Ly=b,U x=y的计算公式；。y =2Z4,？=2,3,：k=lx“=yniunn,七=(/一=-1,一2,i.k=i+l三、程序设计fun c tio n fl=LUre s o lve (a ,b)；n,n =size (a );%行数z=s i z e (b)%b 的行数r；矩阵u=;%u 矩阵fo r j=1 :1 ：nu(l,j)=a(1 ,j);e n dfo r i=2 :1 :n1 (i,l)=a(i ,1)/u(1,1);e n df o r i=2:1:(n 1)sum1=0;fo r k=l:1:(i-1)sum 1 =1 (i,k)*u(k,i)+sum 1 ；e n du(i,i )=a(i,i)-s u ml;s u m 2 =0;s um 3 =0 ；fo r j=(i +1)：1：nfo r k=1 :1 ：(i 1)s u m2=sum2+1 (i,k)*u(k,j )；sum3 =sum 3 +1 (j,k)*u(k,i)e n du(i ,j)=a(i,j)-s u m2;1 (j,i)=(a (j,i)-sum3)/u(i,i);e n de n dfo r i=l:1 ：(n-1 )1 (i,i)=1 ；1 (i,n)=0 ；e n d1 (n,n)=1 ;s u m 4 =0；fo r k=l:1:(n 1)s u m 4=s u m 4+l (n ,k)*u(k,n);e n du(n,n)=a (n,n)s um4;1%输出1矩阵u%输出u矩阵y=l b%输出向量yx=u y%输出向量x a=4 2-3-1 2 1 0 0 0 08 6-5 -3 6 5 0 1 0 04 2-2-1 3 2-1 0 3 10-2 1 5 -1 3-1 1 9 4-4 2 6-1 6 7-3 3 2 38 6-8 5 7 17 2 6-3 50 2-1 3-4 2 5 3 0 116 10-11-9 17 34 2-1 2 24 6 2 -7 1 3 9 2 0 1 2 40 0-1 8-3-24-8 6 3-1 ；b=5;12;3;2;3；46;13；38;1 9；-2 1 ；LUresolve(a,b)z 二10 11 =1.0 0 0 00000000002.000 01.000000000000001.0 0 00001.0 0 000000000-1.0 0 00-2.00 0 03.00 0 01.000 000001.0 0004.0000-0.4 6 671.0000000002.0 0 0 0 1.0000-5.0 0 0 0-0.2667 1.64411.0 0 0 0000 00-1.0000 3.0 0 00-0.06 6 7-0.08 4 7 0.2 9 6 91.00000004.0000-1.000 06.0000-0.26 6 7-3.76272.7303 3.1 863 1.00 0 00 01.0000-4.00 0 0 2 6.000 01.26 6 7 5.00 0 0-0.0182-18.535610.7592 1.0 0 0 000-7.00003 0.00 0 04.60 0 021.7288-1 1.1 1 48-59.930629.2 1 790.98561.0 0000*1.0e+030.00 2 0-0.0030-0.0 0 1 00.00200.001 00.0 02 00.0 0 1 00.00500.010 00.00700.00200.00200.004 000000.003 00.0 02 000000.0 0 100.0 0 20-0.0030-0.00200.0 0 1 0-0.0 0100000.01500.0 1 3 00.03100.0 4 0 00.054 00.063 00.0 6500000-0.003 90.0 1 5 50.03410.07030.0 927 0.1261000000.04170.05180.1 7280.33610.56170 0 0-0.267 20 0 01 669-0.34950 0 00.9 987-1.85620 0 00-0.72 2 90 00 0.0062-0.0 297-0.1 2 1 40 0 00-0.0 840-0.0 0 0000 0 00 0l.Oe+03*0.0 05 00.0 0 20-0.0 0 200.01200.0 1 960.05940.0 0 5 8-0.07 180.84281.849 88.3 7 7 84 1.27 141 0.5 27 830.7 3 80-2 8.0 4 387.35 5 0-14.8 4783.72733.91 2 1-2.5589四、结果分析和讨论LU 分解所得到的结果 x=(8 .3 778,4 1.2 7 1 4 ,1 0.5 2 7 8,3 0.7 3 8 0 ,-2 8.0 4 3 8 ,7.3 5 5 0,-1 4.8 4 7 8,3.72 73,3.9 1 2 1,-2.5 5 8 9 ),与精确解具有彳艮大的误差。这是由于系数矩阵自身的数值性质所决定的，由于计算过程中并未进行选主元的过程，所以由系数矩阵产生的L和U矩阵就具有了很大的数值问题。由L和U所求出的x和y就会产生很大的误差。这是由矩阵自身的数值问题所引起的，与算法自身无关，消除误差的关键在于计算过程中需要进行选主元。五、完毕题目的体会与收获对于其他解线性方程组的方法来看，L U分解是较为高效的，理解并纯熟运用L U分解对于学习数值计算与解线性方程组都有很大的帮助。在进行分解的过程中应注意矩阵的数值问题与如何选取主元的问题，这样才干得到方程组的精确解，否则将产生非常大的误差。在进行分解时应当格外注意，由于系数矩阵存在很多的数值问题。(2)设对称正定阵系数阵线方程组x*=(1,-1,0,2,1-1,0,2/-42-402400%-0-22-1-21320 6-4-1141-8-356了3200-216-1-4-332321-8-1224-10-3九5943-3-44111-4-22025-3-101142入7-150063-3-4219_X8_45 _二、数学原理1、平方根法解n 阶线性方程组Ax=b 的ch o lesk ly 方法也叫做平方根法,这里对系数矩阵A是有规定的，需要A是对称正定矩阵,根据数值分析的相关理论，假如A对称正定,那么系数矩阵就可以被分解为的A=LU 形式，其中L是下三角矩阵,将其代入Ax=b中，可得：LLTx=b进行如下分解：y=Ux x=p f p f(A,b)x=1 2 1.1 4 8 1-1 4 0.1 1 2 72 9.75 1 5-6 0.1 5 2 81 0 .9 1 2 0-2 6.7 9 6 35.4 2 5 9-2.0 1 8 52、改善平方根法fu n c t i o n x=i m p r o v ech o 1 e s k y(A,x=bL =z er o s (n,n)；D=d i a g (n,0);S=L *D;f o r i=l：nL (i,i )=1;en dfo r i =1：nf o r j=l:ni f(e i g (A)n =8;x=i m p r o v e ch o l es k y(A,b,n)x=1 2 1.1 4 8 1-1 4 0.1 1 2 72 9.75 1 5-6 0.1 5 2 81 0.9 1 2 0-2 6 .79 6 35.4 2 5 9-2.0 1 8 5四、结果分析和讨论平方根法和改善平方根法求解线性方程组的解为x=(1 2 1.1 4 8 1,-1 4 0.1 1 2 7,2 9.751 5,-6 0.1 5 2 8,1 0.9 1 2 0,-2 6.7 9 6 3,5 .4 2 5 9,-2.0 1 8 5 )与精确解相比较也存在很大的误差，虽然系数矩阵的对角元素都大于零，原则上可以不必选择主元，但由于矩阵的数值问题较大，不选主元的结果就是产生很大的误差，所以在求解的过程中还是应当选择主元以此消除误差，提高精度。五、完毕题目的体会与收获对称正定矩阵的平方根法及改善平方根法是目前解决这类问题的最有效的方法之一，合理运用的话,可以产生很好的求解效果。改善平方根法较平方根法，由于不用进行开方运算，所以具有一定的求解优势。通过求解此题，使我学会了如何运用Matlab编程来解决平方根法和改善平方根法问题。(3)三对角形线性方程组4-100000000-14-100000000-14-100000000-14-100000000-14-100000000-14-100000000-14-100*X(2,1,-3,0,1,-2,3,0,1,-1/二、数学原理设系数矩阵为三对角矩阵*q0 -000、a2b2C2 一 000A=044,.00000o.an-如%0o 0a.bn 则方程组Ax=f称为三对角方程组。设矩阵A非奇异，A有C rou t分解A=LU,其中L为下三角矩阵,U为单位上三角矩阵,记400 00、1%0.00、a2A0 0001r2.00L=0制A -00,u=00i.00000 Pn-0000 0、000%Be00.1可先依次求出L,U中的元素后，令 U x=y,先求解下三角方程组L y=f得出y,再求解上三角方程组U x =y o事实上，求解三对角方程组的2追赶法将矩阵三角分解的计算与求解两个三角方程组的计算放在一起，使算法更为紧凑。其计算公式为：P P对 i =2,3,%=生 0 i=ba/T，力 唱(*)对 i =,为=%一%为+1三、程序设计f u n c t i o n x =c h a s e(a,b,c ,f)外求解线性方程组A x=f,其中A是三对角阵%a 是矩阵A的下对角线元素a(l)=0%b 是矩阵A的对角线元素枇是矩阵A的上对角线元素c(n)=0%f是方程组的右端向量n=leng th(f);x=z er o s (1 ,n);y =z er o s (1,n);d=z e r o s (1 ,n)；u=z e r o s (1,n);%预解决d (1)=b (1);fo r i=l:n 1u (i)=c(i)/d (i);d (i+l)=b (i+1 )-a(i+1 )*u(i );e n d为追的过程y (1 )=f(D/d(l);fo r i=2 :ny(i )=(f(i)-a(i)*y (i-l)/d (i);en d%赶的过程x (n)=y (n);f or i=n 1 ：-l:1x(i)=y(i)-u(i)*x(i+l);end a=0,-1,-1,-1 ,-l;b =4,4,4,4,4,4,4,4,4,4 ;c=-1 ,-1 ,-1,-1,-1,-1 ,-l,1 ,-1,0 ;f=7,5,-1 3,2,6 ,-1 2,1 4 ,-4,5,-5 ;x =c h a s e(a,b,c,f)x2.0 0 0 01.0 0 0 0-3.0 0 0 00 .0 0 0 01.0 0 0 0-2.0 0 0 03.0 0 0 0-0.0 0 0 01 .0 0 0 0-1.0000四、结果分析和讨论追赶法求解的结果为x=(2,1,-3,0,1,-2,3,0,1,-1)、求解结果与精确解同样，这表白追赶法对于求解三对角方程组具有非常高的精度，误差非常小。算法次数也较少，不选主元也可以有效的算出精确结果,是一种计算量少而数值稳定的方法。五、完毕题目的体会与收获通过本题的求解,加深了对追赶法求解三对角方程组的算法原理的理解。运用追赶法的M a t la b编程,并学会了又一种求解特殊方程组的方法。在计算量方面，追赶法有着巨大的优势，因此在也许的情况下应优先使用追赶法。加深了对数值计算教材知识的理解,收获非常大。题目四:微分方程数值解在传染病的传染理论中，一个比较初等的微分方程可用于预测在任何时刻人群中受传染的数量，只要做了适当的简化。特别的，假定在一个固定的人群中，所有的人具有同样的机会被感染，且一感染就保持这种状态。假设X。）表达在时刻f易受感染的人的数量,y Q）表达感染别人的人数。有理由假设感染别人的人数变化的速率与x 和y的乘积成正比，由于速率既取决于感染别人的人数也取决于那个时刻易感染的人数。假如人口多的足以假定x 和y是连续的的变量，则问题表达为yt）=kx（t）y（t），其中 是常数,必）+y ）=m,m即为人口总数。这个方程就变为仅包含 （/）：问题:一个地区假定m=1 0 0 0 0 0,y（O）=l 0 0 0次=2 x 1 0：又假定期间用天来衡量，求3 0天结束时感染别人的人口数量的近似值二、数学原理E u 1 a r方法：一阶线性微分方程初值问题y=f（x,y）,a x bJ（a）=N oa=x0 xt .川=%,+/矿（X。，%）通过初始值光，依据递推公式（2）逐步算出“以，.,”就为显性的E u l a r方法。隐形E u l a r方法:%=%+/矿(%，M)(3JHI=%+好 区+1，”+1)公式(3)即为隐式E u 1 a r公式。三、程序设计f u n c ti o n E 2=E u l e r_2(f u n,x O,y O,x N,N)%向后E u l e r公式,其中，枇u n为一阶微分方程的函数%x 0,y O为初始条件%x N为取值范围的一个端点助为区间步长额为区间的个数炮 为x N构成的向量%y为yN构成的向量x =z e ro s(1,N+1 );y=z e ro s(l,N +1 );x(l)=x 0 ;y(l)=y 0;h =(x N-x O)/N;f o r n=l:N为用迭代法求y (n+l)x (n+l)=x(n)+h;z O=y (n)+h*f e v a l(f u n,x (n),y (n)；f o r k=1 :3z 1 =y (n)+h*f e v a 1 (f u n,x(n+l),z 0);i f a bs(z 1-z O)l e-3bre a k;e n dz 0=zl；endy(n+l)=zl;en dT=xz,y f u ncti o n z=f(x,y)z=(2*10 (-6)*(1 0000 0-y)*y；Euler_2 C 0,1 0 00,29,29)T=1.0e+04*00.1 00 00.0 0 0 10.12460.0 0020.155 10.00030.19 2 90.00 0 40.239 60.00 0 50.2 97 20.0 0 060.3 6 8 00.000 70.45480.000 80.56 0 50.0 0090.6 8860.0 0100.8 4290.0 01 11.027 10.00 1 20.00130.0 0 1 40.0 0 10.00 1 60.00170.0 0180.00190.0 02 00.0 0 210.00220.002 30.0 0 2 40.00250.0 0 260.00270.00280.0 0 291.2 4 501.49991.79435 2.12 9 52.5 0492.91823.3 6 473.83774.3 28 74.8 2815.32605.81 286.28006.72087.1 3 007.504 57.84288.1 4 49四、结果分析和讨论用隐式欧拉法求得第3 0天时感染别人的人口数量的近似值为y(29)=8 1 44 9人。隐式欧拉法较之欧拉法的误差更小,求取的结果具有一定的可参考价值。对于结果存在的误差，只能通过改变算法实现。五、完毕题目的体会与收获求解初值问题，欧拉法是最简朴的数值解法。而隐式欧拉法更进一步减少了欧拉法的误差。隐式欧拉法的原理依旧十分的简朴。在对数值结果的精确度并非规定很高时可采用隐式欧拉法,这样既能快速地得到结果,又能得到合适的结果。本题的解决对我数值计算的学习有着很大的帮助，使我对初值问题的理解更加深刻，并熟悉了初值问题的求解方法和求解过程,。题目五:数值积分在光学物理中研究矩角位的光绕射经常会用到F r e s n e l 积分：c(t)=c o s(y co2 和 s(f)=s i n(co2)da)对于，=0.1,0.2,1.0构造一个精确到1()7的c(f)和s Q)值的表格供研究者查询(R。m b e rg积分算法)。二 数学原理龙贝格算法是由递推算法得来的。由梯形公式得出辛普生公式得出柯特斯公式最后得到龙贝格公式。在变步长的过程中探讨梯形法的计算规律。设将求积区间 a,b 分 为 n个等分，则一共有n+1 个等分点4=。+山,/=，no这里用7；表达复化梯形法求得的积分值，其n下标n表达等分数。先 考 察 下 一 个 字 段 其 中 点、=；(4+X*M),在该子段上二分前后两个积分值k+2 2I.1=5 1/(/)+)灯，、(、,(=了 f M+f xt i+/(%)4L k k+2)J显然有下列关系+=聂+-将这一关系式关于k从 0 到 n-1累加求和，即可导出下列递推公式,1 ”=尹+工/X，A=0 I 2 7需要强调指出的是，上 式 中 的 心 代 表 二 分 前 的 步 长，而n梯形法的算法简朴，但精度低，收敛速度缓慢，如何提高收敛速度以节省计算量，自然式人们极为关心的。根据梯形法的误差公式,积分值7；的截断误差大体与小成正比，因此步长减半后误差将减至四分之一，既有f 4将上式移项整理,知f产 2-%)由此可见，只要二分前后两个积分值7；和耳相称接近，就可以保证计算保证结果计算结果七的误差很小，这种直接用计算结果来估计误差的方法称作误差的事后估计法。按上式，积 分 值 耳 的 误 差 大 体 等 于 7；),假如用这个误差值作为耳的一种补偿,可以盼望,所得的 1 4 1T-T2 n+-(T2n-Tn)-T2n-T 应当是更好的结果。按上式,组合得到的近似值 直接验证,用梯形二分前后的两个积分值7；和 Q 按式组合，结果得到辛普生法的积分值S,。再考察辛普生法。其截断误差与川成正比。因此，若将步长折半，则误差相应的减至十六分之一。既有f 1 6由此得/S,-S1 5 2 1 5”不难验证，上式右端的值其实就等于C,就是说，用辛普生法二分前后的两个积分值S“和 S*,在按上式再做线性组合,结果得到柯特斯法的积分值C,既有16 1 0C Q 3 c-o15 2 15 反复同样的手续，依据斯科特法的误差公式可进一步导出龙贝格公式A 64 1 。3。3应当注意龙贝格公式已经不属于牛顿一柯特斯公式的范畴。在步长二分的过程中运用公式加工三次,就能将粗糙的积分值7；逐步加工成精度较高的龙贝格凡，或者说,将收敛缓慢的梯形值序列Tn加工成纯熟迅速的龙贝格值序列Rn,这种加速方法称龙贝格算法。三 程序设计f u nc t i on 1 b g(f x,a,t,k,e)%fx为规定的积分函数%a为规定的积分的下限%t 为规定的积分的上限%k 为规定的积分的最大次数%e为规定积分的结束精度T=z e ros (k,k);%T为龙贝格算法递推表T(l,1)=(t-a)*(l +f x(t)/2;f or i =l :km=0;f or j=1:2*(i-1)m=m+f x(a+(2*j T)*(t-a)/(2i);e ndT(i+1,1)=0.5*T(i,l)+(t-a)*m/2 i ;f or n=1:iT(i +1 ,n+1)=(4-n*T (i+1,n)-T (i,n)/(4 n-1);e ndi f a b s(T(i+l ,i+l)-T (i ,i)=4b r e a k；e l s ee n de ndf or i=1 :ki f T(i,1 )=0J =i；b re a k;e l s ee nde ndi f j=ke rror(7所求次数不够或不可积)e 1 s ee ndT=T（1d i s p （所求的积分值为：）%输出求得的积分值d i s p(T(j-1,1)f u n c t i o n f x =f (x )f x=c os(p i *x-2/2)；f u n c t i on f x =f (x )f x=s i n(p i*x八 2/2);a=0 ;t =0.1 ;k=l 0 0;e=1 0(4);f x=(x)co s(pi*x-2/2)；1 b g(f x ,a,t,k,e)T =0.1 0 0 000000.1 0 0 00.1 0 0 000 00.1 0 0 0 0.1 0 0 00.1 0 0 0000.1 0 0 00.1 0 0 00.1 0 0 00.1 0 0 000.1 0 0 00.1 0 0 00.1 0 0 00.1 0 0 00.1 0 0 0所求的积分值为：0.1 0 0 0 a =0；t=0.2;k=1 0 0 ；e =1 0 八(-4);f x=(x )c o s (pi*x -2/2);l b g (f x,a ,t,k,e )T =0 .1 9 9 800000.1 9 9 90.1 9 9 90000.1 9 9 90.1 9 9 90 .1 9 9 9000.1 9 9 90 .1 9 9 90.1 9 9 90.1 9 9 900.1 9 9 90.1 9 9 90.1 9 9 90.1 9 9 90.1 9 9 9所求的积分值为：0.1 9 9 9 a =0 ;t =0.3;k=1 0 0 ;e =1 0 (-4)；f x=(x)c os(pi*x -2/2)；Ib g (f x,a,t,k,e)T =0.2 9 8 500000.2 9 9 20.2 9 9 40000.2 9 9 30.2 9 9 40.29 9 4000.29 9 40.2 9 9 40.29 9 4 0.2 9 9 400.2 9 9 40.29 9 40 .2 9 9 40.2 9 9 40.29 9 4所求的积分值为：0.2 9 9 4 a=0 ;t=0.4;k=1 0 0 ；6=1 0 人(-4)；f x =(x)c os(p i*x -2/2);l b g(f x ,a ,t,k,e)T =0.3 9 37 0 0 0 00.39 6 5 0.39 7 4 0 0 00.3 9 7 2 0.3 9 7 5 0.3 9 7 5 0 00.39 7 4 0.3 9 7 5 0.39 7 5 0.39 7 5 00.39 7 5 0.39 7 5 0.3 9 7 5 0.39 7 5 0.39 7 5所求的积分值为：0.39 7 5 a =0;t =0.5;k=1 0 0;e=1 0 (-4)；f x=(x )c o s (pi*x 人 2/2);l b g(f x,a,t,k,e )0 .57 9 20.48 1 000000.48 9 30.4 9 210000.4 9 1 60.49 2 30.49 2 3000.49 2 10.4 9 2 30.49 2 30.49 2 300.4 9 2 3所求的积分值为：0.49 230.4 92 30.4 9 2 30.49 230.4 9 23 a=0;t=0.6;k=1 0 0;e=1 0 (-4)1 b g(f x,a,t,k,e )T 二;f x=(x)C 0s (p i *x -2/2);0.5 5 3 300000.5 7 370.58 0 40000.58 1 10 .5 8 1 1000.58 0 6 0.58 1 1 0.58 1 1 0.58 1 1 00.58 1 0 0.58 1 1 0.58 1 1 0.5 8 1 1 0.58 1 1所求的积分值为：0.58 1 0 a=0 ；t=0.7;k=1 0 0;e=1 0 7-4)；f x=(x)c os (pi*x 2/2)；1 b g (f x,a,t,k,e)T =0.6 0 1 300000.6 4420 .6 5 8 50000.6 5580.6 59 60.6 59 7000.6 58 70.6 59 60.6 59 70 .6 59 700.6 59 40.6 5 9 70.6 59 70.6 5 9 70.6 59 7所求的积分值为：0.6 5 9 4 a=0；t=0.8 ;k=1 0 0 ;e=l 0 八(-4);f x =(x)c o s (p i*x-2/2);1 b g(f x,a ,t,k,e)T =0.6 1 4300000.6 9 460 .7 2 1 40000.7 1 580.7 2 280.7 229000.7 21 10.7 2280.7 2 2 80.7 22800.7 2 24 0.7 2 2 80.7 2 280.7 2 2 80.7 228所求的积分值为:0.7 224 a=0;t =0.9 ;k=1 0 0;e =1 0*(4);f x=(x)c os (pi /2)；1 b g (f x,a,t,k,e)0.58 2 300000.7 1 8 60.7 6 4 00000.7 5 340 .7 6 5 0 0.7 6 5 0000.7 6 2 00.7 6 4 80.7 6 4 80.7 6 4800.7 6 41 0.7 6 48 0.7 6 4 8 0 .7 6 48 0.7 6 4 8所求的积分值为：0.7 6 4 1 a=0；t=l；k=1 0 0 ；e =1 0 (-4)；f x=(x)c o s (pi*x 2/2);1 b g (f x ,a,t,k,e)T =0.7 7 8 9 0.7 7 9 90.5 00 000000.7 11 90.7 8 2 60000.76 340.7 8 0 50.7 8 0 4000.7 7 580 .7 7 9 90 .7 7 9 90.7 7 9 900.7 7 9 90.7 7 9 90.7 7 9 9所求的积分值为:0.7 789 a=0;t=0.l;k=100;e=10(-4 );fx=(x)sin(p i*x*2/2)；lbg(f x,a,t,k,e)T=0.050 80000000000.0 25 60.017200000 00 00.01300.00 8 90.0083000 00 000.00 6 80.0 04 70.00 4 40.0 0 440000 000.0 0 3 60.00 2 60.00250.00240.002 4000 000.00210.001 60.00 1 50.0 0 1 50.00150.00 1 500000.00130.00100.00 1 00.00100.00100.00 1 00.001000 00.00090.0 0 080.00080.00 0 80.0 0 080.00080.00 0 8 0.0 00 8000.0 0 0 70.00 0 70.00060.0 0 060.0 0 06o.oooe 0.000 60.00060.000600.0 0 0 60.0 0060.0 0060.00 0 60.00