龙贝格算法的matlab实现

实验二数值方法计算积分学号：姓名：指导教师：实验目的1、了解并掌握matlab软件的基本编程、操作方法；2、初步了解matlab中的部分函数，熟悉循环语句的使用；3、通过上机进一步领悟用复合梯形、复合辛普森公式，以及用龙贝格求积方法计算积分的原理。

一、用不同数值方法计算积分 10x ln xdx=-94.(1)取不同的步长h.分别用复合梯形及辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能再被改善？(2)用龙贝格求积计算完成问题（1）。

二、实现实验1、流程图：下图是龙贝格算法框图：2、 算法：(1) 复合梯形公式：Tn=++)()([2b f a f h2∑-=11)](n k xk f ；(2) 复合辛普森公式：Sn=6h[f(a)+f(b)+2∑-=11)](n k xk f +4∑-=+1)2/1(n k x f ]；以上两种算法都是将a-b 之间分成多个小区间（n ）,则h=(b-a)/n,x k =a+kh, x k+1/2=a+(k+1/2)h,利用梯形求积根据两公式便可。

(3) 龙贝格算法：在指定区间内将步长依次二分的过程中运用如下公式1、Sn=34T2n-31Tn 2、 Cn=1516S2n-151Sn3、 Rn=6364C2n-631Cn 从而实现算法。

3、 程序设计（1）、复合梯形法：function t=natrapz(fname,a,b,n) h=(b-a)/n;fa=feval(fname,a);fb=feval(fname,b);f=feval(fname,a+h:h:b-h+0.001*h);t=h*(0.5*(fa+fb)+sum(f));（2）、复合辛普森法：function t=natrapz(fname,a,b,n)h=(b-a)/n;fa=feval(fname,a);fb=feval(fname,b);f1=feval(fname,a+h:h:b-h+0.001*h);f2=feval(fname,a+h/2:h:b-h+0.001*h); t=h/6*(fa+fb+2*sum(f1)+4*sum(f2));（3）龙贝格法：function [I,step]=Roberg(f,a,b,eps)if(nargin==3) eps=1.0e-4;end; M=1; tol=10; k=0;T=zeros(1,1); h=b-a;T(1,1)=(h/2)*(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b));while tol>epsk=k+1;h=h/2;Q=0;for i=1:Mx=a+h*(2*i-1);Q=Q+subs(sym(f),findsym(sym(f)),x);endT(k+1,1)=T(k,1)/2+h*Q;M=2*M;for j=1:kT(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1); endtol=abs(T(k+1,j+1)-T(k,j));endI=T(k+1,k+1);step=k;4、实验结果；（1）复合梯形法(2)、复合辛普森法结果：(3)龙贝格法结果四．总结由结果（1）、（2）可知复合辛普森法求积分精度明显比复合梯形法求积的精度要高，且当步长取不同值时即n 越大、h 越小时，积分精度越高。

matlab中用龙贝格外推法计算积分使用龙贝格外推法计算积分是一种常用的数值积分方法，它通过将积分区间不断细分，并利用数值逼近的思想来求解积分值。

在Matlab中，我们可以使用龙贝格外推法来计算各种类型的积分，从而得到比较准确的结果。

我们需要了解什么是数值积分。

在数学中，积分是函数与自变量之间的一种运算关系，它描述了函数在一定区间上的累积效应。

而数值积分则是通过将积分区间进行离散化，将连续的积分变为离散的求和，从而通过计算得到近似的积分值。

龙贝格外推法是一种常用的数值积分方法，它通过不断细分积分区间，并利用数值逼近的思想来逐步提高积分的精度。

具体来说，龙贝格外推法的计算过程如下：1. 首先，我们将积分区间[a, b]等分为n个小区间，其中n是一个初始值。

2. 然后，我们计算出每个小区间上的积分值。

在Matlab中，可以使用数值积分函数如quad、quadl或quadgk来计算每个小区间的积分值。

3. 接下来，我们使用这些积分值来构造一个n阶的龙贝格外推表。

具体来说，我们可以使用以下公式计算第k+1阶的外推值：R(k+1,1) = 1/(4^k-1) * (4^k * R(k,1) - R(k-1,1))其中，R(k,1)表示第k阶的外推值。

4. 我们重复步骤3，直到达到我们期望的精度要求或迭代次数达到限制。

在每一次迭代中，我们将n加倍，并重新计算每个小区间上的积分值和外推值。

通过上述步骤，我们可以得到一个逐步接近真实积分值的序列。

通常情况下，我们会设定一个误差容限，当两次迭代的结果之差小于该容限时，就认为计算结果已经足够精确。

需要注意的是，龙贝格外推法在计算积分时需要进行一定的迭代计算，因此在计算复杂积分时可能需要较长的计算时间。

此外，为了保证计算结果的准确性，我们还需要选择合适的积分区间和适当的迭代次数。

龙贝格外推法是一种常用的数值积分方法，它通过将积分区间不断细分，并利用数值逼近的思想来求解积分值。

数值分析——数值积分实习题管理科学与工程学院 学号：1120140500 姓名：彭洋洋 一、计算实习题1.用不同数值方法计算积分：049xdx =-⎰.（1）取不同的步长h ，分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h 的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h ，使得精度不能再被改善？ （2）用龙贝格求积计算完成问题（1） （3）用自适应辛普森积分，使其精度达到10-4解答：（1）取不同的步长，采用不同的公式，比较精度过程如下： 1.1 复合梯形公式及复合辛普森公式求解复合梯形公式：11*[()2()()]2n n k k hT f a f x f b -==++∑误差关于h 的函数：2(2)()**()12n a b R f h f ξ-=复合辛普森公式：111/201*[()4()2()()]6n n n k k k k hS f a f x f x f b --+===+++∑∑误差关于h 的函数：4(4)()*(/2)*()180n a bR f h f η-=1.2 复合梯形公式及复合辛普森公式Matlab 程序（2）用龙贝格求积计算完成问题（1） 2.1 龙贝格求积算法龙贝格求积公式也称为逐次分半加速法。

它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。

作为一种外推算法，它在不增加计算量的前提下提高了误差的精度。

24133n n n S T T =- 21611515n n n C S S =- 26416363n n n R C C =-1221/201()22n n n k k h T T f x -+==+∑ ()(1)()11(4*)/(41)k m k k mm m m T T T +--=-- 1,2,...k = 2.2 龙贝格求积Matlab 程序画图程序设计①得到关于n各种公式求积的图表如下：对于梯形公式、辛普森公式n代表份数，龙贝格公式n表示从1开始的序列号②关于步长h 的各种公式求积的图表如下其中龙贝格序列步长()/2k h b a =-：观察两幅图表h 越小，精度越高。

北京科技⼤学应⽤计算⽅法作业与答案⼀、第⼀次作业（⼀）2-6计算下列向量的1-范数、∞-范数、2-范数。

（1）x=(12,-4,-6,2)T >> A=[12,-4,-6,2] A =12 -4 -6 2 >> norm(A,1) ans = 24>> norm(A,inf) ans = 12>> norm(A,2) ans =14.1421 (2) x=(1,3,-4)T >> A=[1,3,-4] A =1 3 -4 >> norm(A,1) ans = 8>> norm(A,inf) ans = 4>> norm(A,2) ans =5.0990（⼆）2-9 计算下列矩阵的⾏范数、列范数、谱范数、F 范数。

（1）--=112111113A >> A=[3,-1,1;1,1,1;2,1,-1]A =3 -1 1 1 1 1 2 1 -1 >> norm(A,1) ans = 6>> norm(A,inf)ans = 5>> norm(A,2) ans =3.7888 >> norm(A,'fro') ans =4.4721 （2）R a a a A ∈?-=,00 >> A=[0,1;-1,0] A =0 1 -1 0 >> norm(A,1) ans = 1>> norm(A,inf) ans = 1>> norm(A,2) ans = 1>> norm(A,'fro') ans =1.4142⼆、第⼆次作业⽤⽜顿迭代法求⽅程0133=--x x 在20=x 附近的根。

要求：给成程序和运⾏结果.1、⽜顿法的基本原理在求解⾮线性⽅程0)(=x f 时，它的困难在于)(x f 是⾮线性函数，为克服这⼀困难，考虑它的线性展开。

多种数值积分求积公式的分析比较吴春晖（中国海洋大学海洋环境学院山东青岛 266100）摘要：对于运用牛顿-莱布尼茨积分公式不能较好解决的定义在区间[a,b]上的可积函数，原函数并不能简单地用初等函数来表达，故需要构造定积分的近似计算公式。

在本文中，主要构建了抛物线求积公式及其复化抛物线公式。

在对抛物线类的求积公式进行应用检验后，再运用Gauss求积公式，构建Gauss-Laguerre求积公式，对相同的问题进行运用，并比较截断误差。

之后再对求积过程进行优化，在限定误差范围的情况下，利用龙贝格算法，对求积加速收敛。

关键词：抛物线求积复化求积 Gauss-Laguerre 加速收敛引言：对于一些较为复杂的函数，在一定的误差要求下，需要通过构造的方式求给定函数的定积分。

基本的替代法主要有梯形面积及抛物线近似代替曲边梯形。

并通过划分更小的区间，减少截断误差从而提出了复化梯形及抛物线公式。

为了提高运算效率，有加速收敛的Richardson外推法和Romberg求积公式。

之后，针对节点数固定情况下，提出了Gauss公式，使其获得最大的精度。

本文主要研究的是抛物线求积法与Gauss-Laguerre公式。

目录第一章抛物线求积公式及应用 (3)1.1抛物线求积公式的算法 (3)1.2抛物线求积公式的matlab程序 (3)1.3复化抛物线求积公式的应用 (4)第二章Gauss-Laguerre求积公式及应用 (5)2.1 Gauss-Laguerre的算法 (5)2.2Gauss-Laguerre公式的matlab程序 (5)2.3 Gauss-Laguerre求积公式的应用 (6)第三章龙贝格算法与算法优化 (7)3.1龙贝格算法及程序 (7)3.2利用龙贝格算法优化求积 (7)3.3 龙贝格算法的应用 (8)第四章数值积分的分析总结 (9)第一章 抛物线求积公式及应用1.1 抛物线求积公式的算法抛物线求积公式，是将区间二等分，以中点及两端点作为抛物线的三个点，并求出抛物线，在区间上对抛物线函数求积分。

《数值分析》课程实验报告一、实验目的1、进一步熟悉向量矩阵的运算；2、掌握龙贝格（Romberg ）算法，并能用高级程序语言MATLAB 编写实现此算法的程序；3、进而加深对龙贝格（Romberg ）算法的理解。

二、实验内容1. 使用Romberg 积分，对于计算下列⎰+4802)cos (1dx x 各近似值a.确定1,51,41,31,21,1,,,,R R R R Rb.确定5,54,43,32,2,,,R R R Rc.6,65,64,63,62,61,6,,,,,R R R R R Rd.确定10,109,98,87,7,,,R R R R三、实验步骤1. 编写程序龙贝格积分方法如下：n=5;a=0;b=48;h(1,1)=b-a;fa=sqrt(1+(cos(a))^2);fb=sqrt(1+(cos(b))^2);r(1,1)=h(1,1)/2*(fa+fb);disp('R11，R21，R31，R41，R51分别为');disp(r(1,1));for i=2:nh(i,1)=(b-a)/(2^(i-1));sum=0;for k=1:2^(i-2)x=a+(2*k-1)*h(i,1);sum=sum+sqrt(1+(cos(x)).^2);endr(i,1)=0.5*(r(i-1,1)+h(i-1,1)*sum);disp(r(i,1));enddisp('R22，R33，R44，R55分别为');for k=2:nfor j=2:kr(k,j)=r(k,j-1)+(r(k,j-1)-r(k-1,j-1))/(4^(j-1)-1);enddisp(r(k,k));enddisp('R61，R62，R63，R64，R65，R66分别为');n=6;for i=2:nh(i,1)=(b-a)/(2^(i-1));sum=0;for k=1:2^(i-2)x=a+(2*k-1)*h(i,1);sum=sum+sqrt(1+(cos(x)).^2);endr(i,1)=0.5*(r(i-1,1)+h(i-1,1)*sum);endfor k=2:nfor j=2:kr(k,j)=r(k,j-1)+(r(k,j-1)-r(k-1,j-1))/(4^(j-1)-1);endendfor i=1:ndisp(r(6,i));enddisp('R77，R88，R99，R10,10分别为');n=10;for i=2:nh(i,1)=(b-a)/(2^(i-1));sum=0;for k=1:2^(i-2)x=a+(2*k-1)*h(i,1);sum=sum+sqrt(1+(cos(x)).^2);endr(i,1)=0.5*(r(i-1,1)+h(i-1,1)*sum);endfor k=2:nfor j=2:kr(k,j)=r(k,j-1)+(r(k,j-1)-r(k-1,j-1))/(4^(j-1)-1);endendfor i=7:10disp(r(i,i));end运行结果如下：R11，R21，R31，R41，R51分别为62.437457.288656.443856.263156.2188R22，R33，R44，R55分别为55.572356.201556.205656.2041R61，R62，R63，R64，R65，R66分别为58.362759.077359.268959.317559.329759.3328R77，R88，R99，R10,10分别为58.422158.470758.470558.4705四、实验小结在这次编程中我学到了很多东西，把程序写进软件中也出现了很多错误，细节问题使我们必须注意的，自己有了很多的收获，自己进一步理解和学习了Matlab软件。

数值分析实验报告实验3 数值积分实验目的通过本实验理解数值积分与微分的基本原理。

掌握数值积分中常见的复合求积公式的编程实现。

掌握龙贝格算法的基本思路和迭代步骤；培养编程与上机调试能力。

算法描述3.2.1 龙贝格算法基本思路（1）将区间[]b a ,划分为n 等分，分点：()n x x 0；根据梯形公式()()()⎥⎦⎤⎢⎣⎡+-=∑-=1122n k k n b f x f a f h T ，求出n T ，再根据n T 和n T 2之间的递推公式()∑-=++=1021222n k k n n x f h T T 求出n T 2； （2）设m 为加速次数，k 为划分区间次数，则由加速公式：()()()k m m k m m m k m T T T 11141144-----=（ ,2,1=k ）求出第k 次划分，第m 次加速次数的梯形值()k m T ，这样不断地循环，直到求出在满足精度条件下的某个()k m T 作为积分值为止。

3.2.2 龙贝格算法计算步骤1.输入MATLAB 程序function[t]=romberg(f,a,b,e)t=zeros(15,4);t(1,1)=(b-a)/2*(f(a)+f(b));for k=2:4sum=0;for i=1:2^(k-2)sum=sum+f(a+(2*i-1)*(b-a)/2^(k-1));endt(k,1)=*t(k-1,1)+(b-a)/2^(k-1)*sum;for i=2:kt(k,i)=(4^(i-1)*t(k,i-1)-t(k-1,i-1))/(4^(i-1)-1); endendfor k=5:15sum=0;for i=1:2^(k-2)sum=sum+f(a+(2*i-1)*(b-a)/2^(k-1));endt(k,1)=*t(k-1,1)+(b-a)/2^(k-1)*sum;for i=2:4t(k,i)=(4^(i-1)*t(k,i-1)-t(k-1,i-1))/(4^(i-1)-1); endif k>6if abs(t(k,4)-t(k-1,4))<edisp(['答案 ',num2str(t(k,4))]);break;endendendif k>=15disp(['溢出']);end2.输入f=inline('sin(x)/x','x')f =Inline function:f(x) = sin(x)/x>> romberg(f,10^(-100),1,5*10^(-7)) 3.运行结果为实验内容用龙贝格算法计算：1 0sin xI dxx=⎰实验步骤3.4.1 代码程序：function[t]=romberg(f,a,b,e)t=zeros(15,4);t(1,1)=(b-a)/2*(f(a)+f(b));for k=2:4sum=0;for i=1:2^(k-2)sum=sum+f(a+(2*i-1)*(b-a)/2^(k-1));endt(k,1)=*t(k-1,1)+(b-a)/2^(k-1)*sum;for i=2:kt(k,i)=(4^(i-1)*t(k,i-1)-t(k-1,i-1))/(4^(i-1)-1); endendfor k=5:15sum=0;for i=1:2^(k-2)sum=sum+f(a+(2*i-1)*(b-a)/2^(k-1));endt(k,1)=*t(k-1,1)+(b-a)/2^(k-1)*sum;for i=2:4t(k,i)=(4^(i-1)*t(k,i-1)-t(k-1,i-1))/(4^(i-1)-1);endif k>6if abs(t(k,4)-t(k-1,4))<edisp(['答案 ',num2str(t(k,4))]);break;endendendif k>=15disp(['溢出']);end3.4.2 实验结果实验体会本次试验使我认识到了计算机计算能力的强大，通过本次实验对数值积分与微分的基本原理有了深刻理解。

桂林理工大学考试（考查）试卷（ 学年度 第一学期）课程名称：数值分析考核专业班级：学号： 姓名： 成绩：﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍﹍一、 简答题（4分×5=20分）1、病态线性方程组的主要判断依据有哪些？答：（1）系数矩阵的某两行（列）几乎近似相关（2）系数矩阵的行列式的值很小（3）用主元消去法解线性方程组时出现小主元（4）近似解x*已使残差向量r=b-Ax*的范数很小，但该近似解仍不符合问题要求。

2、数值计算中如何避免不稳定的算法，防止有效数字的损失？答：（1）简化计算过程，减少运算次数；（2）避免两个相近的数相减；（3）避免除数的绝对值远小于被除数的绝对值；（4）防止大数“吃掉”小数的现象；（5）使用数值稳定的算法，设法控制误差的传播。

3、解释龙格现象。

答：增加插值节点，提高插值多项式的次数，可以使插值函数在更多的点与所逼近的函数取相同的值，但会使插值函数在两端发生激烈的振荡，这就是插值计算的龙格现象。

4、阐述迭代法的基本思想。

答：就是用某种极限过程去逐步逼近线性方程组精确解得方法。

其基本思想为：先任取一组近似解初值X 0，然后按照某种迭代原则，由X 0计算新的近似解X 1，以此类推，可计算出X 2,X 3，…X K ，。

,如果｛X ｝收敛，则取为原方程组的解。

5、叙述QR 算法的步骤。

答：对于给定的n 阶实对称矩阵A 与迭代次数M ；（1）令A 1=A ，对于k=1,2，…，M ；（2）迭代计算下一个矩阵：A k =Q k R k （对A k 作QR 分解），（3）A k =RQ(交换乘法次序)，令k=k+1,A k+1=Q k T A K Q K（4）返回到2，直到k=M,输出A 的主对角元素。

二、 计算型编程题（20分×3=60分）1、 用Newton 插值拟合[0, 2*pi]上函数()sin x f x e x 。

答：function [yt,N] = NewtInterp(x,y,xt)% 已知数据点的牛顿插值% x,y:插值条件% xt:要计算的插值点，可以是多个% yt:用牛顿插值函数出xt对应的函数值数组% N: 牛顿插值多项式表达式syms t;n=length(x);ny=length(y);if n~=nyerror('插值节点x与函数值y维数不一致'); enda=zeros(1,n);N = y(1);w = 1;for k=1:n-1yy=zeros(1,n); % 记录差商for j=k+1:nyy(j) = (y(j)-y(k))/(x(j)-x(k));enda(k) = yy(k+1);w = w*(t-x(k));N = N + a(k)*w;y = yy;endyt = subs(N,'t',xt);simplify(N);N = collect(N); % 将插值多项式展开N = vpa(N, 6); % 系数转化为6位精度输入：x=[0:pi/10:2*pi];y=exp(x).*sin(x);xt=pi;[yt,N] = NewtInterp(x,y,xt)% 画图z=0:pi/20:2*pi;yz= subs(N,'t',z); %计算插值点处的函数值figure;plot(z,exp(z).*sin(z),'--r',z,yz,'-b')hold onplot(x,y,'marker','+')hold onplot(xt,yt,'marker','o')h=legend('$\exp{x}.*sin{x}$','Newton','$(x_k,y_k)$','$x=pi$');set(h,'Interpreter','latex')xlabel('x')ylabel('y')yt =7.1684e-015N =.159313e-14*t^20-.178421e-12*t^19+.578696e-11*t^18-.993977e-10*t^17+.1179 00e-8*t^16-.101872e-7*t^15+.640214e-7*t^14-.340920e-6*t^13+.131789e-5*t^12-.33 7754e-5*t^11+.19392e-4*t^10+.22800e-4*t^9+.33873e-4*t^8-.162956e-2*t^7-.110707 e-1*t^6-.333620e-1*t^5+.1611e-4*t^4+.333328*t^3+1.00000*t^2+.999992*t2、利用龙贝格算法程序计算积分x。