数值计算B大作业Word版

1、课 程 设 计课程名称： 数值计算B 设计题目： 数值计算B大作业 学 号： 姓 名： 完成时间： 整理为word格式题目一：多项式插值某气象观测站在8：00（AM）开始每隔10分钟对天气作如下观测，用三次多项式插值函数（Newton）逼近如下曲线，插值节点数据如上表，并求出9点30分该地区的温度（x=10）。x12345678y22.523.324.421.7025.228.524.825.4二、数学原理假设有n+1个不同的节点及函数在节点上的值（x，y），（x，y），插值多项式有如下形式： （1）其中系数（i=0,1,2n）为特定系数，可由插值样条（i=0,1,2n）确定。根据均差的定义，

2、把x看成a,b上的一点,可得f(x)= f（）+f（）fx, = f+fx, （）fx, ，x= fx, ，x+ fx, ，x(x-x)综合以上式子，把后一式代入前一式，可得到： f(x)= f+f（）+ f（）（）+ fx, ，x（）(x-x)+ fx, ，x,= N（x）+其中N（x）= f+f（）+ f（）（）+ fx, ，x（）(x-x) （2）= f(x)- N（x）= fx, ，x, （3） 整理为word格式=（）(x-x)Newton插值的系数（i=0,1,2n）可以用差商表示。一般有 （k=0,1,2，n ） （4）把（4）代入（1）得到满足插值条件N（i=0,1,2，n）的

3、n次Newton插值多项式N（x）=f（）+f（）+f（）（）+f（）（）（）.其中插值余项为： 介于之间。三、程序设计function y,A,C,L=newdscg(X,Y,x,M) % y为对应x的值，A为差商表，C为多项式系数，L为多项式% X为给定节点，Y为节点值，x为待求节点n=length(X); m=length(x); % n为X的长度for t=1:m z=x(t); A=zeros(n,n);A(:,1)=Y's=0.0; p=1.0; q1=1.0; c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)- A(i-1,j-1)

4、/(X(i)-X(i-j+1); end q1=abs(q1*(z-X(j-1);c1=c1*j; end C=A(n,n);q1=abs(q1*(z-X(n);for k=(n-1):-1:1C=conv(C,poly(X(k); d=length(C);C(d)=C(d)+A(k,k);end y(k)= polyval(C, z); %输出y值endL(k,:)=poly2sym(C); %输出多项式整理为word格式>> syms M,X=1,3,5,7;Y=22.5,24.4,25.2,24.8;x=10;>> y,A,C,L=newdscg(X,Y,x,M)

5、y = 21.7313A = 22.5000 0 0 0 24.4000 0.9500 0 0 25.2000 0.4000 -0.1375 0 24.8000 -0.2000 -0.1500 -0.0021C = -0.0021 -0.1187 1.4521 21.1688 L = - x3/480 - (19*x2)/160 + (697*x)/480 + 3387/1604、 结果分析和讨论 对于不超过三次的插值多项式，x如果选取1,3,5,7这三个点能够得到较好的三次插值多项式L=-0.0021x3-0.1187x2+1.4521x+21.1688。当x=10时，也即9点30分时的温度

6、为21.7317度，结果分析知此值应是偏小的。对于选取不同的插值节点，能够得到不同的插值多项式，误差也不尽相同。5、 完成题目的体会与收获 牛顿插值法的重要一点就是对插值节点的选取，通过本题的编程很好的加深了对其概念的理解以及提高了应用牛顿插值法的能力，学会了运用Matlab软件对牛顿插值法相关问题进行编程求解，对Matlab计算方法与程序编辑更加熟悉。使我对这类问题的理解有了很大的提升整理为word格式。整理为word格式题目二：曲线拟合在某钢铁厂冶炼过程中，根据统计数据的含碳量与时间关系，试用最小二乘法拟合含碳量与时间t的拟合曲线，并绘制曲线拟合图形。t(分)0 5 10 15 20 25

7、 30 35 40 45 50 55 0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64 2、 数学原理从整体上考虑近似函数同所给数据点(i=0,1,m)误差(i=0,1,m)的大小，常用的方法有以下三种：一是误差(i=0,1,m)绝对值的最大值，即误差 向量的范数；二是误差绝对值的和，即误差向量r的1范数；三是误差平方和的算术平方根，即误差向量r的2范数；前两种方法简单、自然，但不便于微分运算 ，后一种方法相当于考虑 2范数的平方，因此在曲线拟合中常采用误差平方和来 度量误差(i=0，1，m)的整体大小。数据拟合的具体作法是：对给

8、定数据 (i=0,1,，m)，在取定的函数类中,求,使误差(i=0,1,m)的平方和最小，即 =从几何意义上讲，就是寻求与给定点(i=0,1,m)的距离平方和为最小的曲线。函数称为拟合 函数或最小二乘解，求拟合函数的方法称为曲线拟合的最小二乘法。在曲线拟合中，函数类可有不同的选取方法。三、程序设计>> x=0,5,10,15,20,25,30,35,40,45,50,55;>> y=0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.02,4.64*10(-4);>> p=polyfit(x,y,2)p =整理为

9、word格式 1.0e-04 * -0.0024 0.2037 0.2305>> plot(x,y,'r')4、 结果分析和讨论 通过最小二乘法得到的曲线拟合多项式是p=(-0.0024x2+0.2037x+0.2305)*10-4。利用Matlab软件调用最小二乘法函数即可得到多项式拟合函数，由于数据较少得到的拟合曲线不够光滑，可能会存在一定的误差。拟合曲线总体呈现增函数趋势，与数据较为吻合。5、 完成题目的体会与收获 曲线拟合较常用的方法之一就包括最小二乘法，因此能够熟练使用Matlab进行数据的曲线拟合变得尤为重要。通过完成本题的拟合过程，对于最小二乘法曲线拟

10、合的操作更加的熟练，能够较好的完成各类数据的拟合。题目三：线性方程组求解整理为word格式分别利用LU分解；平方根法与改进平方根法；追赶法求解下列几个不同类型的线性方程组，并与准确值比较结果，分析误差产生原因。（1）设线性方程组 二、数学原理 将A分解为一个下三角矩阵L和一个上三角矩阵U，即：A=LU,其中 L=， U= 解Ax=b 的问题就等价于要求解两个三角形方程组： Ly=b,求y;   Ux=y,求x。 设A为非奇异矩阵,且有分解式A=LU, L为单位下三角阵,U为上三角阵。 L,U的元素可以有n步直接计算定出。用直接三角分解法解Ax=b（要求A的所有顺序主子式都不

11、为零）的计算公式： , ,i=2,3,，n. 计算U的第r行，L的第r列元素(i=2,3,n)：     ， i=r,r+1,n; , i=r+1,n,且rn.整理为word格式 求解Ly=b，Ux=y的计算公式； 三、程序设计function f1=LUresolve(a,b);n,n=size(a); % 行数z=size(b) % b的行数l=; % l矩阵u=; % u矩阵for j=1:1:n u(1,j)=a(1,j);endfor i=2:1:n l(i,1)=a(i,1)/u(1,1);endfor i=2:1:(n-1) sum1=0; for

12、k=1:1:(i-1) sum1=l(i,k)*u(k,i)+sum1; end u(i,i)=a(i,i)-sum1; sum2=0; sum3=0; for j=(i+1):1:n for k=1:1:(i-1) sum2=sum2+l(i,k)*u(k,j); sum3=sum3+l(j,k)*u(k,i); end u(i,j)=a(i,j)-sum2; l(j,i)=(a(j,i)-sum3)/u(i,i); endendfor i=1:1:(n-1) l(i,i)=1; l(i,n)=0;endl(n,n)=1;整理为word格式sum4=0;for k=1:1:(n-1) sum

13、4=sum4+l(n,k)*u(k,n);endu(n,n)=a(n,n)-sum4;l %输出l矩阵u %输出u矩阵y=lb %输出向量yx=uy %输出向量x>> a=4 2 -3 -1 2 1 0 0 0 0 8 6 -5 -3 6 5 0 1 0 0 4 2 -2 -1 3 2 -1 0 3 1 0 -2 1 5 -1 3 -1 1 9 4 -4 2 6 -1 6 7 -3 3 2 3 8 6 -8 5 7 17 2 6 -3 5 0 2 -1 3 -4 2 5 3 0 1 16 10 -11 -9 17 34 2 -1 2 2 4 6 2 -7 13 9 2 0 12 4

14、 0 0 -1 8 -3 -24 -8 6 3 -1;>> b=5;12;3;2;3;46;13;38;19;-21;>> LUresolve(a,b)z = 10 1l = 1.0000 0 0 0 0 0 0 0 0 0 2.0000 1.0000 0 0 0 0 0 0 0 0 1.0000 0 1.0000 0 0 0 0 0 0 0 0 -2.0000 3.0000 1.0000 0 0 0 0 0 0 -1.0000 1.0000 4.0000 -0.4667 1.0000 0 0 0 0 0 2.0000 1.0000 -5.0000 -0.2667 -1

15、.6441 1.0000 0 0 0 0 0 -1.0000 3.0000 -0.0667 -0.0847 0.2969 1.0000 0 0 0 4.0000 -1.0000 6.0000 -0.2667 -3.7627 2.7303 3.1863 1.0000 0 0 1.0000 -4.0000 26.0000 1.2667 5.0000 -0.0182 -18.5356 10.7592 1.0000 0 0 -7.0000 30.0000 4.6000 21.7288 -11.1148 -59.9306 29.2179 0.9856 1.0000u =整理为word格式 1.0e+03

16、 * 0.0040 0.0020 -0.0030 -0.0010 0.0020 0.0010 0 0 0 0 0 0.0020 0.0010 0.0050 0.0100 0.0070 0.0020 0.0030 0.0020 0.0020 0 0 0.0010 0 0.0020 0 -0.0030 -0.0020 0.0010 -0.0010 0 0 0 0.0150 0.0130 0.0310 0.0400 0.0540 0.0630 0.0650 0 0 0 0 -0.0039 0.0155 0.0341 0.0703 0.0927 0.1261 0 0 0 0 0 0.0417 0.05

17、18 0.1728 0.3361 0.5617 0 0 0 0 0 0 0.0062 -0.0297 -0.1214 -0.2672 0 0 0 0 0 0 0 -0.0840 -0.1669 -0.3495 0 0 0 0 0 0 0 0 -0.9987 -1.8562 0 0 0 0 0 0 0 0 0 -0.7229y = 1.0e+03 * 0.0050 0.0020 -0.0020 0.0120 0.0196 0.0594 0.0058 -0.0718 0.8428 1.8498x = 8.3778 41.2714 10.5278 30.7380 -28.0438 7.3550 -1

18、4.8478 3.7273 3.9121 -2.5589整理为word格式4、 结果分析和讨论 LU分解所得到的结果x=(8.3778,41.2714,10.5278,30.7380,-28.0438,7.3550,-14.8478,3.7273,3.9121,-2.5589)T与精确解具有很大的误差。这是由于系数矩阵本身的数值性质所决定的，由于计算过程中并未进行选主元的过程，所以由系数矩阵产生的L和U矩阵就具有了很大的数值问题。由L和U所求出的x和y就会产生很大的误差。这是由矩阵本身的数值问题所引起的，与算法本身无关，消除误差的关键在于计算过程中需要进行选主元。5、 完成题目的体会与收获对于

19、其他解线性方程组的方法来看，LU分解是较为高效的，理解并熟练运用LU分解对于学习数值计算与解线性方程组都有很大的帮助。在进行分解的过程中应注意矩阵的数值问题与如何选取主元的问题，这样才能得到方程组的精确解，否则将产生非常大的误差。在进行分解时应该格外注意，因为系数矩阵存在很多的数值问题。（2）设对称正定阵系数阵线方程组 2、 数学原理1、 平方根法解n阶线性方程组Ax=b的choleskly方法也叫做平方根法，这里对系数矩阵A是有要求的，需要A是对称正定矩阵，根据数值分析的相关理论，如果A对称正定，那么系数矩阵就可以被分解为的形式，其中L是下三角矩阵，将其代入Ax=b中，可得：进行如下分解：整

20、理为word格式那么就可先计算y,再计算x，由于L是下三角矩阵，是上三角矩阵，这样的计算比直接使用A计算简便，同时你应该也发现了工作量就转移到了矩阵的分解上面，那么对于对称正定矩阵A进行Cholesky分解，我再描述一下过程吧：如果你对原理很清楚那么这一段可以直接跳过的。设，即其中第1步，由矩阵乘法，故求得一般的，设矩阵L的前k-1列元素已经求出第k步，由矩阵乘法得于是2、 改进平方根法在平方根的基础上，为了避免开方运算,所以用计算；其中，；得整理为word格式 按行计算的元素及对元素公式 对于. 计算出的第行元素后，存放在的第行相置，然后再计算的第行元素，存放在的第行.的对角元素存放在的相应

21、位置. 对称正定矩阵按分解和按分解计算量差不多，但分解不需要开放计算。求解, 的计算公式分别如下公式。 3、 程序设计1、平方根法function x=pfpf(A,b)%楚列斯基分解 求解正定矩阵的线性代数方程 A=LL 先求LY=b 再用LX=Y 即可以求出解Xn,n=size(A);L(1,1)=sqrt(A(1,1);for k=2:n L(k,1)=A(k,1)/L(1,1);end整理为word格式for k=2:n-1 L(k,k)=sqrt(A(k,k)-sum(L(k,1:k-1).2); for i=k+1:n L(i,k)=(A(i,k)-sum(L(i,1:k-1).*

22、L(k,1:k-1)/L(k,k); endendL(n,n)=sqrt(A(n,n)-sum(L(n,1:n-1).2);%解下三角方程组Ly=b 相应的递推公式如下，求出y矩阵y=zeros(n,1);%先生成方程组的因变量的位置，给定y的初始值for k=1:n j=1:k-1; y(k)=(b(k)-L(k,j)*y(j)/L(k,k);end%解上三角方程组 LX=Y 递推公式如下，可求出X矩阵x=zeros(n,1);U=L'%求上对角矩阵for k=n:-1:1 j=k+1:n; x(k)=(y(k)-U(k,j)*x(j)/U(k,k);end >> A=4

23、,2,-4,0,2,4,0,0 2,2,-1,-2,1,3,2,0 -4,-1,14,1,-8,-3,5,6 0,-2,1,6,-1,-4,-3,3 2,1,-8,-1,22,4,-10,-3 4,3,-3,-4,4,11,1,-4 0,2,5,-3,-10,1,14,2 0,0,6,3,-3,-4,2,19;>> b=0;-6;20;23;9;-22;-15;45;>> x=pfpf(A,b)x = 121.1481 -140.1127整理为word格式 29.7515 -60.1528 10.9120 -26.7963 5.4259 -2.01852、改进平方根法f

24、unction x=improvecholesky(A,b,n) %用改进平方根法求解Ax=bL=zeros(n,n); %L为n*n矩阵D=diag(n,0); %D为n*n的主对角矩阵S=L*D;for i=1:n %L的主对角元素均为1 L(i,i)=1;endfor i=1:n for j=1:n %验证A是否为对称正定矩阵if (eig(A)<=0)|(A(i,j)=A(j,i) %A的特征值小于0或A非对称时，输出wrong disp('wrong');break;endendendD(1,1)=A(1,1); %将A分解使得A=LDLTfor i=2:n f

25、or j=1:i-1 S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1)'); L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1); end D(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-1)');endy=zeros(n,1); % x，y为n*1阶矩阵x=zeros(n,1);for i=1:n y(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)/D(i,i); %通过 LDy=b解得y的值endfor i=n:-1:1 x(i)=y(i)-sum(

26、L(i+1:n,i)'*x(i+1:n); %通过LTx=y解得x的值end整理为word格式>> A=4,2,-4,0,2,4,0,0 2,2,-1,-2,1,3,2,0 -4,-1,14,1,-8,-3,5,6 0,-2,1,6,-1,-4,-3,3 2,1,-8,-1,22,4,-10,-3 4,3,-3,-4,4,11,1,-4 0,2,5,-3,-10,1,14,2 0,0,6,3,-3,-4,2,19;>> b=0;-6;20;23;9;-22;-15;45;>> n=8;>> x=improvecholesky(A,b,n)

27、x = 121.1481 -140.1127 29.7515 -60.1528 10.9120 -26.7963 5.4259 -2.01854、 结果分析和讨论 平方根法和改进平方根法求解线性方程组的解为x=（121.1481，-140.1127，29.7515，-60.1528，10.9120，-26.7963，5.4259，-2.0185）T。与精确解相比较也存在很大的误差，虽然系数矩阵的对角元素都大于零，原则上可以不必选择主元，但由于矩阵的数值问题较大，不选主元的结果就是产生很大的误差，所以在求解的过程中还是应该选择主元以此消除误差，提高精度。 5、 完成题目的体会与收获 对称正定矩阵

28、的平方根法及改进平方根法是目前解决这类问题的最有效的方法之一，合理利用的话，能够产生很好的求解效果。改进平方根法较平方根法，因为不用进行开方运算，所以具有一定的求解优势。通过求解此题，使我学会了如何运用整理为word格式Matlab编程来解决平方根法和改进平方根法问题。（3）三对角形线性方程组 2、 数学原理 设系数矩阵为三对角矩阵则方程组Ax=f称为三对角方程组。 设矩阵A非奇异，A有Crout分解A=LU，其中L为下三角矩阵，U为单位上三角矩阵，记整理为word格式可先依次求出L，U中的元素后，令Ux=y，先求解下三角方程组Ly=f得出y，再求解上三角方程组Ux=y。 事实上，求解三对角方

29、程组的2追赶法将矩阵三角分解的计算与求解两个三角方程组的计算放在一起，使算法更为紧凑。其计算公式为：(*)三、程序设计function x=chase(a,b,c,f)%求解线性方程组Ax=f,其中A是三对角阵%a是矩阵A的下对角线元素a(1)=0%b是矩阵A的对角线元素%c是矩阵A的上对角线元素c(n)=0%f是方程组的右端向量n=length(f);x=zeros(1,n);y=zeros(1,n);d=zeros(1,n);u= zeros(1,n);%预处理d(1)=b(1);for i=1:n-1u(i)=c(i)/d(i);d(i+1)=b(i+1)-a(i+1)*u(i);end

30、%追的过程y(1)=f(1)/d(1);for i=2:n y(i)=(f(i)-a(i)*y(i-1)/d(i);end%赶的过程x(n)=y(n);整理为word格式for i=n-1:-1:1x(i)=y(i)-u(i)*x(i+1);end>> a=0,-1,-1,-1,-1,-1,-1,-1,-1,-1;>> b=4,4,4,4,4,4,4,4,4,4;>> c=-1,-1,-1,-1,-1,-1,-1,-1,-1,0;>> f=7,5,-13,2,6,-12,14,-4,5,-5;>> x=chase(a,b,c,f)x

31、=2.00001.0000 -3.00000.00001.0000 -2.00003.0000 -0.00001.0000 -1.00004、 结果分析和讨论 追赶法求解的结果为x=（2,1，-3,0,1，-2,3,0,1，-1）T。求解结果与精确解一样，这表明追赶法对于求解三对角方程组具有非常高的精度，误差非常小。算法次数也较少，不选主元也可以有效的算出精确结果，是一种计算量少而数值稳定的方法。5、 完成题目的体会与收获 通过本题的求解，加深了对追赶法求解三对角方程组的算法原理的理解。运用追赶法的Matlab编程，并学会了又一种求解特殊方程组的方法。在计算量方面，追赶法有着巨大的优势，因此在

32、可能的情况下应优先使用追赶法。加深了对数值计算教材知识的理解，收获非常大。题目四：微分方程数值解整理为word格式在传染病的传染理论中，一个比较初等的微分方程可用于预测在任何时刻人群中受传染的数量，只要做了适当的简化。特别的，假定在一个固定的人群中，所有的人具有同样的机会被感染，且一感染就保持这种状态。假设表示在时刻易受感染的人的数量，表示感染别人的人数。有理由假设感染别人的人数变化的速率与和的乘积成正比，因为速率既取决于感染别人的人数也取决于那个时刻易感染的人数。如果人口多的足以假定和是连续的的变量，则问题表示为，其中是常数，即为人口总数。这个方程就变为仅包含：问题：一个地区假定，又假定时间

33、用天来衡量，求30天结束时感染别人的人口数量的近似值2、 数学原理Eular方法： 一阶线性微分方程初值问题 （1）方程离散化：差分和差商 （2） 通过初始值，依据递推公式（2）逐步算出就为显性的Eular方法。隐形Eular方法：整理为word格式 （3）公式（3）即为隐式Eular公式。3、 程序设计function E2=Euler_2(fun,x0,y0,xN,N)% 向后Euler公式，其中，%fun为一阶微分方程的函数%x0,y0为初始条件%xN为取值范围的一个端点%h为区间步长%N为区间的个数%x为xN构成的向量%y为yN构成的向量x=zeros(1,N+1);y=zeros(1

34、,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;for n=1:N%用迭代法求y(n+1) x(n+1)=x(n)+h; z0=y(n)+h*feval(fun,x(n),y(n); for k=1:3 z1=y(n)+h*feval(fun,x(n+1),z0); if abs(z1-z0)<1e-3 break; end z0=z1; end y(n+1)=z1;end T=x',y'function z=f(x,y)z=(2*10(-6)*(100000-y)*y;>> Euler_2('f',0,1000,29,29

35、)T =整理为word格式 1.0e+04 * 0 0.1000 0.0001 0.1246 0.0002 0.1551 0.0003 0.1929 0.0004 0.2396 0.0005 0.2972 0.0006 0.3680 0.0007 0.4548 0.0008 0.5605 0.0009 0.6886 0.0010 0.8429 0.0011 1.0271 0.0012 1.2450 0.0013 1.4999 0.0014 1.7943 0.0015 2.1295 0.0016 2.5049 0.0017 2.9182 0.0018 3.3647 0.0019 3.8377 0

36、.0020 4.3287 0.0021 4.8281 0.0022 5.3260 0.0023 5.8128 0.0024 6.2800 0.0025 6.7208 0.0026 7.1300整理为word格式 0.0027 7.5045 0.0028 7.8428 0.0029 8.14494、 结果分析和讨论 用隐式欧拉法求得第30天时感染别人的人口数量的近似值为y（29）=81449人。隐式欧拉法较之欧拉法的误差更小，求取的结果具有一定的可参考价值。对于结果存在的误差，只能通过改变算法实现。5、 完成题目的体会与收获 求解初值问题，欧拉法是最简单的数值解法。而隐式欧拉法更进一步减少了欧拉

37、法的误差。隐式欧拉法的原理依旧十分的简单。在对数值结果的精确度并非要求很高时可采用隐式欧拉法，这样既能快速地得到结果，又能得到合适的结果。本题的解决对我数值计算的学习有着很大的帮助，使我对初值问题的理解更加深刻，并熟悉了初值问题的求解方法和求解过程，。题目五：数值积分整理为word格式在光学物理中研究矩角位的光绕射经常会用到Fresnel积分：对于构造一个精确到的值的表格供研究者查询（Romberg积分算法）。2、 数学原理 龙贝格算法是由递推算法得来的。由梯形公式得出辛普生公式得出柯特斯公式最后得到龙贝格公式。在变步长的过程中探讨梯形法的计算规律。设将求积区间a，b分为n个等分，则一共有n+

38、1个等分点，n。这里用表示复化梯形法求得的积分值，其下标n表示等分数。先考察下一个字段，其中点，在该子段上二分前后两个积分值显然有下列关系将这一关系式关于k从0到n-1累加求和，即可导出下列递推公式需要强调指出的是，上式中的代表二分前的步长，而梯形法的算法简单，但精度低，收敛速度缓慢，如何提高收敛速度以节省计算量，自然式人们极为关心的。根据梯形法的误差公式，积分值的截断误差大致与成正比，因此步长减半后误差将减至四分之一，既有整理为word格式将上式移项整理，知由此可见，只要二分前后两个积分值和相当接近，就可以保证计算保证结果计算结果的误差很小，这种直接用计算结果来估计误差的方法称作误差的事后估

39、计法。按上式，积分值的误差大致等于，如果用这个误差值作为的一种补偿，可以期望，所得的应当是更好的结果。按上式，组合得到的近似值直接验证，用梯形二分前后的两个积分值和按式组合，结果得到辛普生法的积分值。再考察辛普生法。其截断误差与成正比。因此，若将步长折半，则误差相应的减至十六分之一。既有由此得不难验证，上式右端的值其实就等于，就是说，用辛普生法二分前后的两个积分值和，在按上式再做线性组合，结果得到柯特斯法的积分值，既有重复同样的手续，依据斯科特法的误差公式可进一步导出龙贝格公式应当注意龙贝格公式已经不属于牛顿柯特斯公式的范畴。在步长二分的过程中运用公式加工三次，就能将粗糙的积分值逐步加工成精度

40、较高的龙贝格，或者说，将收敛缓慢的梯形值序列加工成熟练迅速的龙贝格值序列，这种加速方法称龙贝格算法。整理为word格式3、 程序设计function lbg(fx,a,t,k,e)% fx为要求的积分函数% a为要求的积分的下限% t为要求的积分的上限% k为要求的积分的最大次数% e为要求积分的结束精度T=zeros(k,k); % T为龙贝格算法递推表T(1,1)=(t-a)*(1+fx(t)/2;for i=1:k m=0; for j=1:2(i-1) m=m+fx(a+(2*j-1)*(t-a)/(2i); end T(i+1,1)=0.5*T(i,1)+(t-a)*m/2i; fo

41、r n=1:i T(i+1,n+1)=(4n*T(i+1,n)-T(i,n)/(4n-1); end if abs(T(i+1,i+1)-T(i,i)<=e & i>=4 break; else ; endendfor i=1:k if T(i,1)=0 j=i; break; else ; endendif j=k error('所求次数不够或不可积')else ;endT=T(1:j-1,1:j-1)disp('所求的积分值为：') % 输出求得的积分值整理为word格式disp(T(j-1,1)function fx = f(x)fx=

42、cos(pi*x2/2);function fx = f(x)fx=sin(pi*x2/2);>> a=0;t=0.1;k=100;e=10(-4);fx=(x)cos(pi*x2/2);>> lbg(fx,a,t,k,e)T = 0.1000 0 0 0 0 0.1000 0.1000 0 0 0 0.1000 0.1000 0.1000 0 0 0.1000 0.1000 0.1000 0.1000 0 0.1000 0.1000 0.1000 0.1000 0.1000所求的积分值为： 0.1000>> a=0;t=0.2;k=100;e=10(-4)

43、;fx=(x)cos(pi*x2/2);>> lbg(fx,a,t,k,e)T = 0.1998 0 0 0 0 0.1999 0.1999 0 0 0 0.1999 0.1999 0.1999 0 0 0.1999 0.1999 0.1999 0.1999 0 0.1999 0.1999 0.1999 0.1999 0.1999所求的积分值为： 0.1999>> a=0;t=0.3;k=100;e=10(-4);fx=(x)cos(pi*x2/2);>> lbg(fx,a,t,k,e)T = 0.2985 0 0 0 0整理为word格式 0.2992 0.2994 0 0 0 0.2993 0.2994 0.2994 0 0 0.2994 0.2994 0.2994 0.2994 0 0.2994 0.2