龙贝格算法的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 越小时，积分精度越高。

1．准备初值，计算
且 （ 为等份次数）
2．按梯形公式的递推关系，计算
3．按龙贝格公式计算加速值
4．精度控制。对给定的精度 ，若
则终止计算，并取 作为所求结果；否则 ，重复2~4步，直到满足精度为止。
问题
（1）程序运行如下：
I = Romberginterg(inline('x.^2.*exp(x)'),0,1,25,1e-6)
function I = GaussInterg(fun, type, a, b, tol)
% GaussInterg用Gauss型求积公式求积分,具体形式由使用者选取
%
% Synopsis: I = GaussInterg(fun, type, a, b)
% I = GaussInterg(fun, type, a, b, tol)
if nargin < 4
npanel = 25;
end
if nargin < 5
tol = 5e-9;
end
if nargin < 6
flag = 0;
end
T(1,1) = TrapezoidInteg(fun, a, b, npanel); %T0(h) = T(h)
err = 1; %初始化误差值
if nargin < 4
npanel = 25;
end
nnode = npanel + 1; %节点数=段数+ 1
h = (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) );

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 越小，精度越高。

3 / 13
一、任务综述于分析
1.1 任务
1、 根据所掌握龙贝格算法，独立编写程序，实现具体问题的求解等； 2、 自己设计一个多项式，根据设定的代数多项式进行测试，调式程序； 3、 对所设计的程序解决实际问题并分析，撰写分析报告
1.2 分析
1、 首先先明确什么是龙贝格算法，了解他的基本原理，画出流程图，编写程序代码。 2、 自己设计一个多项式，对他利用龙贝格算法求出他的积分。并与真实值比较，检测程序是否有错误。 并进一步修改完善程序。可以多试几个多项式，来判断程序的正确性。 3、 任务三实际上是求平面曲线长的问题。利用弧长微分公式可求得结果。
四、测试 ................................................................................................................................................. 8
主要参考资料
《数值分析简明教程》 《数值分析课程设计》
王能超 陈越、童若锋
高等教育出版社 浙江大学出版社
教研室签字： 年 月 日
院签字： 年 月 日
2 / 13
目录
目录 ......................................................................................................................................................... 3 一、任务综述于分析 ............................................................................................................................. 4

北京科技⼤学应⽤计算⽅法作业与答案⼀、第⼀次作业（⼀）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。

用ＭＡＴＬＡＢ计算某些区域上的二重积分摘要：本文研究某些区域上二重积分数值积分公式的构造及用数学软件MATLAB 实现所构造的数值积分的计算，通过MATLAB 用所得公式计算某些典型的二重积分。

主要工作包括：将定积分数值计算的几个公式推广到二重积分，编制MATLAB 程序，最后通过具体的数值算例进行精度比较，从中选出精度高的二重积分的计算公式，并用公式计算一些典型的二重积分。

关键词：二重积分；数值计算；插值多项式；求积公式 1.引言二重积分的计算在科学计算中起着重要的作用，关于矩形区域上的二重积分的计算一般都是化重积分为累次积分，然后借助定积分已有的数值积分计算公式推导出，MATLAB 已经有这些计算公式的相应的函数，但是往往我们建模得到的二重积分的积分区域都不是矩形区域，对于一般的非矩形区域的二重积分，直接用MATLAB 是无法计算的。

又当被积函数比较复杂，无法用初等函数表示或求其原函数很困难时，就只能求积分的数值解。

若),(y x f 在D 上连续，二重积分y x y x f Dd d ),(⎰⎰存在且为一确定的常数，这个数值与),(y x f 的结构、D 的几何形状有关，二重积分计算的基本途径是在一定条件下化为二次积分，本文研究的某些区域的二重积分，要求二重积分在该区域上能化为二次积分。

二重积分的存在性[1]：),(y x f 在闭区域D 上连续，则(,)Df x y dxdy ⎰⎰必存在。

定理[1]：若),(y x f 在闭区域D {;()()}≤≤≤≤a x b c x y d x 上连续，且()c x 、()d x 在[,]a b 上连续，则：()()(,)(,)b d x ac x DI f x y dxdy f x y dxdy ==⎰⎰⎰⎰上式右端是一个先对x 后对y 的二次积分：先把),(y x f 看作x 的函数，在区间[(),()]c x d x 上对x 计算定积分（这时y 看作常数），把得到的结果（是y 的函数）再在[,]a b 上对y 计算定积分即为二重积分。

页脚内容1实验四 数值积分与数值微分专业班级：信计131班 姓名：段雨博 学号：2013014907 一、实验目的1、熟悉matlab 编程。

2、学习数值积分程序设计算法。

3、通过上机进一步领悟用复合梯形、复合辛普森公式，以及用龙贝格求积方法计算积分的原理。

二、实验题目 P1371、用不同数值方法计算积分049xdx =-⎰。

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

三、实验原理与理论基础1.1复合梯形公式及其复合辛普森求解[]()()()11101()()222n n n k k k k k h h T f x f x f a f x f b --+==⎡⎤=+=++⎢⎥⎣⎦∑∑误差关于h 的函数：()()212n b a R fh f η-''=-页脚内容2复合辛普森公式：()()()()111/201426n n n k k k k h S f a f x f x f b --+==⎡⎤=+++⎢⎥⎣⎦∑∑误差关于h 的函数：()()441802n n b a h R f I S f η-⎛⎫=-=- ⎪⎝⎭1.2龙贝格求积算法：龙贝格求积公式是梯形法的递推化，也称为逐次分半加速法，它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种计算积分的方法，同时它有在不断增加计算量的前提下提高误差的精度的特点。

计算过程如下：（1）取0,k h b a ==-，求：()()()[]()00.,.2hT f a f b k a b =+→⎡⎤⎣⎦令k 1记为区间的二分次数 （2）求梯形值02k b a T -⎛⎫⎪⎝⎭即按递推公式12102122n n n k k h T T f x -+=⎛⎫=+ ⎪⎝⎭∑计算0k T .（3）求加速值，按公式()()()111444141m m k k k mm m m m T T T +--=---逐个求出T 表的地k 行其余各元素()()1,2,,k j j T j k -=（4）若()()001k k T T ε--<（预先给定的精度），则终止计算，并取()()0;1k T I k k ≈+→否则令转（2）继续计算。

MATLAB复化梯形法及龙贝格法计算定积分复化梯形法是一种数值积分方法，用于计算定积分的近似值。

该方法的基本思想是将积分区间等分成多个子区间，并在每个子区间上使用梯形公式来进行近似计算。

具体步骤如下：1.将积分区间[a,b]等分成n个子区间，每个子区间的长度为h=(b-a)/n。

2.在每个子区间上，使用梯形公式计算近似积分值。

梯形公式可以表示为：T=(f(x0)+f(x1))*h/2，其中x0和x1分别是子区间的左右边界，f(x)是被积函数。

3.对所有子区间的近似积分值进行求和，得到整个积分区间的近似积分值。

复化梯形法的精度可以通过增加子区间的数量来提高，即使n越来越大，积分值的近似精度也会越来越高。

以下是一个用MATLAB实现复化梯形法计算定积分的示例代码：```matlabh=(b-a)/nresult = 0;for i = 0:n-1x0=a+i*h;x1=a+(i+1)*h;result = result + (f(x0) + f(x1)) * h / 2;endend```接下来，我们来介绍龙贝格法，龙贝格法是一种迭代数值积分方法，用于计算定积分的近似值。

该方法的基本思想是在梯形公式的基础上应用Richardson外推技术，通过逐步加密和外推，提高积分值的精度。

具体步骤如下：1.初始化一个矩阵，矩阵的第一列为复化梯形法的近似积分值。

2.逐列递推计算，每一列的元素为由前一列的元素计算得到。

计算公式为：R(j,k+1)=R(j,k)+(R(j,k)-R(j-1,k))/((4^k)-1)其中，R(j,k)是第j次迭代中计算的近似积分值，k表示第k次迭代。

3.判断是否达到预设的精度要求，如果满足要求，则返回最终近似积分值；否则，继续迭代计算。

以下是一个用MATLAB实现龙贝格法计算定积分的示例代码：```matlabfunction result = romberg(f, a, b, epsilon, max_iter)R = zeros(max_iter, max_iter);h=b-a;R(1,1)=h*(f(a)+f(b))/2;for k = 2:max_iterh=h/2;sum = 0;for i = 1:2^(k-2)x=a+(2*i-1)*h;sum = sum + f(x);endR(k, 1) = R(k-1, 1) / 2 + h * sum;for j = 2:kR(k,j)=R(k,j-1)+(R(k,j-1)-R(k-1,j-1))/((4^(j-1))-1); endif abs(R(k, k) - R(k-1, k-1)) < epsilonresult = R(k, k);return;endendresult = R(max_iter, max_iter);end```这个代码定义了一个名为`romberg`的函数，它接受五个参数：被积函数`f`、积分区间的左边界`a`、积分区间的右边界`b`、精度要求`epsilon`和最大迭代次数`max_iter`。

各种数值积分方法的MATLAB 程序和分析 --以一个实际问题为例问题回顾从某条河的横截面上获得与河岸不同距离(y,m)处的深度数据(H,m)如下，确定这条河的平问题分析由于所得数据为离散点，并且数据点的间隔不均匀，因此先对离散数据进行连续化处理，再对所得函数进行数值积分。

离散数据点的分布图如下：一般来说，不考虑泥沙沉降因素的河底，用变分法来求取切力最小的河床面，可以知道椭圆是最理想的截面形状，然而根据上图可知，河床底部明显存在泥沙沉降。

并且该河段可能处于流向曲率较大的河段，导致内河床(x <5)泥沙沉降较多，河床较浅；而外河床(x >5)部分由于水流速度大于泥沙启动速度，河床底部泥沙被水流携带，因而较深。

根据泥沙启动概率与流速的关系可将泥沙沉降量与河岸距离近似处理为三次关系，如下图：河岸距离 m河道深度 m泥沙沉降量与离河岸距离的关系离河岸距离泥沙沉降量（负值为携带量)由于泥沙启动速度为止，那么将泥沙沉降和携带量达到平衡的地方作为拟合参数处理。

在理想河床截面形状的基础上以上述三次关系修正形状，则回归的函数模型如下：H=α1√1−α2(y−5)2−α3y(y−α4)(y−10)上式中前一项为理想河床截面的椭圆方程，后一项为泥沙沉降量的修正项，先利用非线性回归对上述数据进行拟合。

非线性回归采用Levenberg–Marquardt算法，其函数fitnonpoly.m 源代码如下：%本函数使用Levenberg–Marquardt算法计算待求数据点在要求函数类下的最佳拟合函数%输入：% x为待拟合数据点的横坐标，y为待拟合数据点的纵坐标% fun为自变量和待求参数的函数句柄，格式为fun(x,a)，其中a为参数向量% n为待求参数向量个数，即a长度% tol为结果要求精度%输出：% funchd为所得的非线性拟合函数，可直接使用function funchd=fitnonpoly(x,y,fun,n,tol)m=length(x);%输入检查if m~=length(y)error('The length of x must equal that of y !')elseif m<=n+1error('The length of x must be larger than n+1 !');end%变量创建a=zeros(1,n)+0.1;Z=zeros(m,n);delt=tol/10; %数值微分的扰动项delta=Inf*ones(1,n); %迭代的参数增量初始化v=5; %阻尼项衰减比lambda=10^-2; %阻尼项初始化flag=0;while norm(delta,Inf)>tol||flag==0 %迭代停止条件为参数增量小于要求精度%计算迭代矩阵方程for p=1:nTemp=zeros(1,n);Temp(p)=delt;Z(:,p)=(fun(x',a+Temp)-fun(x',a))./delt;endD=y'-fun(x',a);%计算迭代所得参数增量向量delta=(Z'*Z+lambda*eye(n))\(Z'*D);if norm(y'-fun(x',a+delta'))<norm(D)%若当前迭代收敛，则更新参数向量a=a+delta';lambda=lambda/v; %衰减阻尼项flag=1; elselambda=lambda*v; %若当前迭代不收敛，则增大阻尼 flag=0; end end%将函数表达式转换为可直接使用的函数句柄 syms tfunchd=@(t)fun(t,a); %依次显示所求得的参数值 disp('The parameter in turn is'); disp(a); end运行附件中脚本文件scriofint.m 可以得到拟合后的函数，拟合结果如下：非线性回归所得函数如下：H =1.782√1−0.3995(y −5)2−0.02014y(y −2.735)(y −10)可见泥沙沉降和携带量达到动态平衡的地方里河岸2.735m ,下面对河岸截面积进行数值积分求取。

目录设计总说明 (1)关键词 (1)绪论 (2)第1章课程设计正文 (3)1.1 需求分析 (3)1.2 界面设计 (3)1.3 功能实现 (4)1.4系统算法和程序代码的实现 (5)第2章总结 (7)参考文献 (8)设计总说明学习了《数值分析》和相关开发软件课程并能够熟练掌握各种算法的原理及其求解方法和特点；1、独立进行面向对象软件的分析、设计及在计算机中进行编写代码和调试工作。

培养自己独立进行面向对象软件的分析、面向对象软件的设计及面向对象软件的实现与测试的能力；2、通过课题的理论设计和在计算机中实验调试代码，加深对数值分析知识的理解，培养软件开发的实践技能，提高分析解决具体问题的能力；3、充分利用面向对象程序设计，用精要的语句设计程序，培养独立创新能力；4、熟练运用MATLAB7开发界面美观、清晰，简单，好学，好用的计算软件；说明书中一般包括任务的提出、方案论证、设计与计算说明、试验调试及结果的分析、结束语等内容。

要求理论依据充分，数据准确，公式推导及计算结果正确。

涉及到计算机软件：学生要独立完成一个软件或较大软件中的一个模块，要有足够的工作量；要写出相关的设计文件；能够进行计算机演示和给出运行结果。

在此软件中，我们运用龙贝格方法进行积分。

此外，只需要在编辑框中输入上下限和含x的被积函数，点击计算按钮，得到计算结果。

关键词：MATLAB 数值分析龙贝格积分绪论数值分析是计算数学的一个主要部分，它研究用计算机求解各种数学问题的数值计算方法及其理论与软件实现。

数值分析一个很强的特点就是面向计算机，要根据计算机特点提供切实可行的有效算法。

因此，本次课程设计是在学习了《数值分析》和相关开发软件课程之后，为了能掌握数值分析的基本概念，结合实际的操作和设计，巩固课堂教学内容，使自己能掌握数值分析的基本概念、原理和技术，将理论与实际相结合，应用现有的数据建模工具和数值分析理论，规范、科学地完成一个小型可视化软件的设计与实现，把理论课与实验课所学内容做一综合，通过实际项目的设计、开发，培养学生独立进行计算软件的建模、在计算机中进行计算、设计、通过相关软件开发系统，并在此基础上强化学生的实践意识、提高其实际动手能力和创新能力。

计算机辅助求解一元高次方程计算机辅助求解一元高次方程高次方程是代数学中的一种常见类型，如二次方程、三次方程、四次方程等，其中最高次数通常不超过四次。

一元高次方程是指只包含一个未知数，并且最高次数为n的代数方程式。

由于高次方程的求解过程比较繁琐，直接使用手算难以完成，因此需要计算机辅助来进行求解，大大提高了求解效率。

本文将介绍一元高次方程的求解过程，以及计算机辅助求解高次方程的方法。

一、一元高次方程的求解过程一元高次方程一般写为：$a_nx^n+a_{n-1}x^{n-1}+…+a_1x+a_0=0$，其中$a_n,a_{n-1},…,a_1,a_0,n$均为已知常数，$n\geqslant 2$。

求解一元高次方程的基本思路是：先通过代数变形将高次项的系数化为1，然后使用某些特定的求根方法求解方程。

1. 代数变形(1) 高次项系数化为1有时为求解方便，可以通过约定高次项系数为1来化简原方程，即：$x^n+a_{n-1}x^{n-1}+…+a_1x+a_0=0$。

这时我们把此方程称为首项系数为1的一元高次方程。

(2) 特殊公式变形对于特殊的一元高次方程，可以通过特殊公式的变形来简化求解过程。

例如：$x^2-2x+1=0$，可以变形为$(x-1)^2=0$，从而得到$x=1$。

$x^2+2x+1=0$，可以变形为$(x+1)^2=0$，从而得到$x=-1$。

(x+1)(x+2)=0，可以得到$x=-1~或~x=-2$。

(3) 用辅助方程消去高次项一元高次方程的求解还可以通过使用辅助方程进行变形求解。

例如对于二次方程：$ax^2+bx+c=0$，可以通过配方法将其变形为：$(x+\frac{b}{2a})^2-\frac{\Delta}{4a^2}=0$，其中$\Delta=b^2-4ac$。

同样的方法，对于三次方程、四次方程等，也可以通过辅助方程简化求解。

2. 求根方法(1) 利用求根公式进行求解对于二次方程，可以通过求根公式：$x_{1,2}=\frac{-b\pm\sqrt{b^2-4ac}}{2a}$进行求解。

龙贝格算法matlab程序一、龙贝格算法简介龙贝格算法是一种数值积分的方法，它可以用来计算函数在给定区间上的定积分。

该算法基于复合梯形公式和复合辛普森公式，通过逐步逼近真实值来得到数值解。

它是一种自适应方法，即在每个子区间上都采用不同的步长以获得更高的精度。

二、龙贝格算法的原理1. 复合梯形公式复合梯形公式是将一个区间分成若干个小区间，在每个小区间上应用梯形公式求出积分值，最后将所有小区间的积分值相加得到整个区间上的积分值。

具体公式如下：$$\int_{a}^{b}f(x)dx \approx \frac{h}{2}(f(a)+2\sum_{i=1}^{n-1}f(a+ih)+f(b))$$其中，h为步长，n为子区间数。

2. 复合辛普森公式复合辛普森公式是将一个区间分成若干个小区间，在每个小区间上应用辛普森公式求出积分值，最后将所有小区间的积分值相加得到整个区间上的积分值。

具体公式如下：$$\int_{a}^{b}f(x)dx \approx \frac{h}{6}(f(a)+4\sum_{i=1}^{n/2-1}f(a+(2i)h)+2\sum_{i=1}^{n/2}f(a+(2i-1)h)+f(b))$$其中，h为步长，n为子区间数。

3. 龙贝格算法龙贝格算法是通过不断加密网格来逼近真实值的方法。

首先，将整个区间分成若干个小区间，并在每个小区间上应用复合梯形公式求出初始积分值T(0,0)，然后将这些积分值相加得到整个区间上的积分值T(0,1)。

接着，在每个小区间上应用复合辛普森公式求出更精确的积分值T(1,0)，并将所有小区间的积分值相加得到整个区间上的积分值T(1,1)。

然后，计算两次结果之差E(1,0)=T(1,1)-T(0,1)，如果E(1,0)小于给定误差限，则直接输出T(1,1)作为积分结果；否则，在每个子区间中再次应用复合辛普森公式求出更精确的积分值，并计算两次结果之差，直到满足误差限为止。