龙贝格积分法C语言(西安交大)

龙贝格公式就是逐次对分积分区间的方法，可以把前面计算的结果作为一个整体带入对分后的计算公式中，只需要增加新的分点的函数值。

以)(0k T表示],[b a 二分k 次后的梯形值，)(k mT 表示)(0k T 的m 次加速值，则有),2,1(,141144)(14)1(1)( =---=-+-k TTTk m k m mmk m龙贝格算法：步1 初始化：计算)]()([2)0(0b f a f ab T+-=，置1:=k （k记录区间],[b a 的二分次数），a b h -=；步2（二分） 计算积分值： ∑-=+---+=122/11)1(0)(01)(221k j j k k k x f h TT；步3（加速） 求加速值：计算)(1)1(1)(141144j k j jj k j jjj k jTTT--+------=，),,2,1(k j =；步4 精度检验：对指定的精度ε，若ε<--)0(1)0(k kTT，则终止计算，并取)0(kT作为所求的结果；否则置1:+=k k ，hh 21:=，转步2。

计算次序：(0)0T第1次循环 二分：(1)0T加速：(0)1T 第2次循环二分：(2)1T加速：(1)(0)12,T T第3次循环 二分：(3)2T加速：(2)(1)(0)123,,T T T第4次循环 二分：(4)3T加速：(3)(2)(1)(0)1234,,,T T TT……在命令窗口输入： a = 0; b = 1;epsilon = 5e-6;f = @(x)sin(x);%@(X)申请变量空间,计算sin 积分y = romberg(f,a,b,epsilon) %函数调用安回车，出结果。

计算方法（C）目录第1章绪论1。

1 数值计算1。

2 数值方法的分析1.2.1计算机上数的运算1.2.2算法分析第2章线性代数方程组2。

1 Gauss消去法2.1.1消去法2.1.2主元消去法2.2 矩阵分解2.2.1Gauss消去法的矩阵意义2.2.2矩阵的LU分解及其应用2.2.3其他类型矩阵的分解2.2.4解三对角矩阵的追赶法2.3线性方程组解的可靠性2.3.1向量和矩阵范数2.3.2残向量与误差的代数表征2.4解线性方程组解的迭代法2.4.1基本迭代法2.4.2迭代法的矩阵表示2.4.3收敛性第3章数据近似3。

1 多项式插值3.1.1插值多项式3.1.2Lagrange插值多项式3.1.3Newton插值多项式3.1.4带导数条件的插值多项式3.1.5插值公式的余项3. 2 最小二乘近似3.2.1 最小二乘问题的法方程3.2.2 正交化算法第4章数值微积分4.1 内插求积，Newton-Cotes公式4.1.1Newton—Cotes公式4.1.2复化求积公式4.1.3步长的选取4.1.4Romberg方法4.1.5待定系数法4.2数值微分4.2.1插值公式方法4.2.2Taylor公式方法 (待定系数法)4.2.3外推法第5章非线性方程求解5。

1 解一元方程的迭代法5.1.1简单迭代法5.1.2Newton法5.1.3割线法5.1.4区间方法5。

2 收敛性问题5.2.1简单迭代—-不动点5.2.2收敛性的改善5.2.3Newton法的收敛性5.2.4收敛速度第1章绪论1。

1数值计算现代科学的发展，已导致科学与技术的研究从定性前进到定量，尤其是现代数字计算机的出现及迅速发展，为复杂数学问题的定量研究与解决，提供了强有力的基础.通常我们面对的理论与技术问题，绝大多数都可以从其物理模型中抽象出数学模型，因此,求解这些数学模型已成为我们面临的重要任务.一、本课程的任务：寻求解决各种数学问题的数值方法—-如何将高等数学的问题回归到初等数学（算术)的方法求解—-了解计算的基础方法，基本结构(否则只须知道数值软件）——并研究其性质.立足点:面向数学——解决数学问题面向计算机—-利用计算机作为工具充分发挥计算机的功能，设计算法，解决数学问题例如：迭代法、并行算法二、问题的类型1、离散问题:例如,求解线性方程组bAx= -—从离散数据:矩阵A和向量b，求解离散数据x；2、连续问题的离散化处理:例如,数值积分、数值微分、微分方程数值解；3、离散问题的连续化处理：例如，数据近似,统计分析计算；1.2数值方法的分析在本章中我们不具体讨论算法，首先讨论算法分析的基础--误差.一般来讲，误差主要有两类、三种（对科学计算）：1）公式误差-—“截断误差”,数学↔计算，算法形成——主观(人为）：数学问题-数值方法的转换,用离散公式近似连续的数学函数进行计算时，一般都会发生误差，通常称之为“截断误差”；——以后讨论2)舍入误差及输出入误差——计算机,算法执行—-客观（机器）：由于计算机的存储器、运算器的字长有限，在运算和存储中必然会发生最末若干位数字的舍入，形成舍入误差；在人机数据交换过程中，十进制数和二进制数的转换也会导致误差发生,这就是输入误差.这两种误差主要是由于计算机的字长有限,采用浮点数系所致。

第1章 龙贝格算法研究问题阐述计算二重积分⎰⎰=b adcdxdy y x f I ),(------------------------------------(1-1)对工科生要求计算如下3个二重积分：算例1：1110()1I x y dxdy =+=⎰⎰算例2：/2200()I x y dxdy πππ=+=⎰⎰-----------------------(1-2)算例3： 1130sin()x y I dxdy x y+=+⎰⎰通常将积分区间等分数依次取2k ，并得到相应的递推式子如（2-5）所示，相应对到过程不再累述详情可以参考教材。

1112221[()()],21[(21)]222k k k k k i b a T f a f b b a b a T T f a i --=-⎧=+⎪⎪⎨--⎪=++-⎪⎩∑----------------------------(2-5) 最终得到计算机计算步骤如下：(1)计算初值T1； (2)K 1；(3)计算新的T2，一直迭代计算2kT ;(4)精度控制：如果122||kk T T ε--<，则停止计算，并将2k T 作为积分的近似数值，否则循环将不停止。

对第四点，通常可以加入一个循环最高次数，防止程序陷入死循环。

function f=fxy(x,y,w) if w==1f=x+y; %函数1elseif w==2f=sin(abs(x-y)); %函数2elseif w==3f=sin(x+y)/(x+y); %函数3if x+y==0f=1;end;end;function gy=gy(y,a,b,w,error)k=1;lock=0;count=0;T1=(b-a)*(fxy(a,y,w)+fxy(b,y,w))/2;while lock==0&&count<100 %限定最多迭代20次T2=T2kx(y,T1,a,b,k,w);if abs(T2-T1)<errorgy=T2;lock=1;elseT1=T2;count=count+1;k=k+1;end;end% count%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 龙贝格算法%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%% 作者：XXXX %%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%% 时间：2012/4/16 %%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%clear all; clc;%% 详情见报告Romberg公式介绍部分clear all;clc;which=1; %表示计算的函数f1,f2,f3%积分区间a=0;if which==1||which==3b1=1;b2=1; %函数1和2的区间elseb1=pi;b2=pi/2; %函数2的区间end%% 计算主程序k=1;lock=0;count=0; %辅助量error=0.0001; %计算精度% gy=gy(y,a,b,w,error)T1=(b2-a)*(gy(a,a,b1,which,error)+gy(b2,a,b1,which,error))/2; while lock==0&&count<100 %限定最多迭代20次T2=T2ky(T1,a,b1,b2,k,which,error);if abs(T2-T1)<errorI=T2;lock=1;elseT1=T2;function T2kx=T2kx(y,T2kx_1,a,b,k,w)temp1=0;for i=1:2^(k-1)temp1=temp1+fxy(a+(2*i-1)*(b-a)/2^k,y,w);endT2kx=0.5*T2kx_1+temp1*(b-a)/2^k;function T2ky=T2ky(T2ky_1,a,b1,b2,k,w,error)temp1=0;for i=1:2^(k-1)temp1=temp1+gy(a+(2*i-1)*(b2-a)/2^k,a,b1,w,error);endT2ky=0.5*T2ky_1+temp1*(b2-a)/2^k;。

龙贝格积分1. 算法原理采用复化求积公式计算时，为使截断误差不超过ε，需要估计被积函数高阶导数的最大值，从而确定把积分区间[]b a ,分成等长子区间的个数n 。

首先在整个区间[]b a ,上应用梯形公式，算出积分近似值T1；然后将[]b a ,分半，对 应用复化梯形公式算出T2；再将每个小区间分半，一般地，每次总是在前一次的基础上再将小区间分半，然后利用递推公式进行计算，直至相邻两个值之差小于允许误差为止。

实际计算中，常用ε≤-n n T T 2作为判别计算终止的条件。

若满足，则取n T f I 2][≈；否则将区间再分半进行计算，知道满足精度要求为止。

又经过推导可知，∑=-++=ni i i n n x x f h T T 112)2(221，在实际计算中，取kn 2=，则k a b h 2-=，112)1*2(2++--+=+k i i ab i a x x 。

所以，上式可以写为∑=++--+-+=+kk i k k ab i a f a b T T 211122)2)12((2211k开始计算时，取())()(21b f a f ab T +-=龙贝格算法是由递推算法得来的。

由梯形公式得出辛普森公式得出柯特斯公式最后得到龙贝格公式。

根据梯形法的误差公式，积分值n T 的截断误差大致与2h 成正比，因此步长减半后误差将减至四分之一，即有21114n n T T -≈-将上式移项整理，知2211()3n n n T T T -≈-由此可见，只要二分前后两个积分值n T 和2n T 相当接近，就可以保证计算保证结果计算结果2n T 的误差很小，这种直接用计算结果来估计误差的方法称作误差的事后估计法。

按上式，积分值2n T 的误差大致等于21()3n n T T -，如果用这个误差值作为2n T 的一种补偿，可以期望，所得的()222141333n n n n n T T T T T T =+-=-应当是更好的结果。

龙贝格（Romberg ）求积法1.算法理论Romberg 求积方法是以复化梯形公式为基础，应用Richardson 外推法导出的数值求积方法。

由复化梯形公式 )]()(2)([2222b f h a f a f h T +++=可以化为)]()]()([2[212112h a f h b f a f hT +++=＝)]([21211h a f h T ++一般地，把区间[a,b]逐次分半k －1次，（k ＝1,2,……,n ）区间长度（步长）为kk m a b h -=，其中mk ＝2k －1。

记k T ＝)1(k T由)1(k T ＝]))12(([21211)1(1∑=---++km j k k k h j a f h T 从而⎰badxx f )(＝)1(kT－)(''122k f h a b ξ- （1）按Richardson 外推思想，可将（1）看成关于k h ，误差为)(2k h O 的一个近似公式，因而，复化梯形公式的误差公式为⎰badxx f )(－)1(k T ＝......4221++kkh K h K ＝∑∞=12i i k i h K （2）取1+k h ＝k h 21有⎰badxx f )(－)1(1+k T＝∑∞=+121221i i k iihK （3）误差为)(2jh O 的误差公式 )(j kT＝)1(-j kT＋141)1(1)1(------j j k j k T T 2.误差及收敛性分析（1）误差，对复化梯形公式误差估计时，是估计出每个子区间上的误差，然后将n 个子区间上的误差相加作为整个积分区间上的误差。

（2）收敛性，记h x i =∆，由于∑=++=ni i i n x f x f h f T 01))]()([2)(＝))()((21101∑∑-==∆+∆n i ni i i i i x x f x x f上面两个累加式都是积分和，由于)(x f 在区间],[b a 上可积可知，只要],[b a 的分划的最大子区间的长度0→λ时，也即∞→n 时，它们的极限都等于积分值)(f I 。

实验三龙贝格方法【实验内容】1.理解龙贝格方法的基本思路2.用龙贝格方法设计算法，编程求解一个数值积分的问题。

【实验方法】实验程序#include<stdio.h>#include<math.h>#define E 0.0000005 //函数近似值的精度float f(float x) //定义f(x)函数{float y;y=4/(1+x*x); //f(x),运行前输入return y;}main(){printf("被求积分的函数为4/(1+X*X) \n \n");float a,b,h;float T1,T2,S1,S2,C1,C2,R1,R2;printf("请分别输入积分的下限a 和积分的上限b ，用空格分开\n");scanf("%f %f",&a,&b); //输入函数上下限h=b-a;T1=0.5*h*(f(a)+f(b)); //T1int k,i;for(i=2;1;i++) //循环计算T S C R{R1=R2;float s=0; //∑F(X(i+0.5))for(k=1;k<=pow(2,(i-2));k++){s=s+f(a+(k-0.5)*h);}T2=0.5*T1+0.5*h*s;S2=(4*T2-T1)/3;T1=T2; //计算完S2后，把T2赋值给T1C2=(16*S2-S1)/15;S1=S2; //计算完C2后，把S2赋值给S1R2=(64*C2-C1)/63;C1=C2; //计算完R2后，把C2赋值给C1h=h/2;if((fabs(R2-R1))<E) //达到精度要求时，退出循环{break;}}printf("函数的近似值为%f \n",R2); //输出}【实验结果】。

龙贝格算法及应用龙贝格算法是一种数值计算方法，用于计算数值积分的近似解。

它的应用范围广泛，涉及到物理、工程、金融等领域。

在下面的回答中，我将详细介绍龙贝格算法的原理、步骤和一些应用实例。

原理：龙贝格算法是一种基于复化梯形公式的数值积分算法。

它通过不断提高插值点数目，逐步提高数值积分的精度。

步骤：1. 首先，我们选择一个足够小的初始步长h，然后根据复化梯形公式计算积分的近似值I(h)。

2. 接下来，我们将步长h减半，再次计算积分的近似值I(h/2)。

3. 然后，根据复化梯形公式，通过I(h/2)和I(h)计算出一个更高阶的数值积分近似值I(2h)。

4. 重复上述步骤，每次将步长减半，直到达到所需的精度。

应用：龙贝格算法广泛应用于数值积分问题，特别是对于某些复杂函数，无法通过解析方法求得精确解的情况下。

以下是一些具体的应用实例：1. 物理学中的轨道运动描述：龙贝格算法可以用来计算行星围绕太阳的轨道运动，以及其他天体运动的数值积分问题。

2. 工程学中的电路分析：通过对电路中的电流和电压的积分，可以计算电路中的功率、能量等物理量。

3. 金融学中的期权定价：龙贝格算法可以用于计算期权定价模型中的积分，以估计期权的市场价值。

4. 数理统计学中的概率密度函数估计：通过计算概率密度函数的积分，可以对数据的分布进行建模和估计。

5. 计算机图形学中的曲线绘制：龙贝格算法可以用来计算曲线的长度、面积、甚至绘制曲线的插值。

总结：龙贝格算法是一种常用的数值积分算法，通过逐步提高插值点数目，能够对于复杂函数进行积分近似。

它在物理、工程、金融等领域有广泛的应用。

通过应用龙贝格算法，我们可以获得数值计算问题的近似解，从而解决实际问题中的积分计算需求。

龙贝格（Romberg ）求积法1.算法理论Romberg 求积方法是以复化梯形公式为基础，应用Richardson 外推法导出的数值求积方法。

由复化梯形公式 )]()(2)([2222b f h a f a f h T +++=可以化为)]()]()([2[212112h a f h b f a f hT +++=＝)]([21211h a f h T ++一般地，把区间[a,b]逐次分半k －1次，（k ＝1,2,……,n ）区间长度（步长）为kk m a b h -=，其中mk ＝2k －1。

记k T ＝)1(k T 由)1(k T ＝]))12(([21211)1(1∑=---++km j k k k h j a f h T 从而⎰badxx f )(＝)1(kT－)(''122k f h a b ξ- （1）按Richardson 外推思想，可将（1）看成关于k h ，误差为)(2k h O 的一个近似公式，因而，复化梯形公式的误差公式为⎰badxx f )(－)1(k T ＝......4221++k k h K h K ＝∑∞=12i i k i h K （2）取1+k h ＝k h 21有 ⎰ba dx x f )(－)1(1+k T ＝∑∞=+121221i ik ii hK （3）误差为)(2jh O 的误差公式 )(j kT＝)1(-j kT＋141)1(1)1(------j j k j k T T2.误差及收敛性分析（1）误差，对复化梯形公式误差估计时，是估计出每个子区间上的误差，然后将n 个子区间上的误差相加作为整个积分区间上的误差。

（2）收敛性，记h x i =∆，由于∑=++=ni i i n x f x f h f T 01))]()([2)(＝))()((21101∑∑-==∆+∆n i ni i i i i x x f x x f上面两个累加式都是积分和，由于)(x f 在区间],[b a 上可积可知，只要],[b a 的分划的最大子区间的长度0→λ时，也即∞→n 时，它们的极限都等于积分值)(f I 。

龙贝格算法 一、问题分析1.1龙贝格积分题目要求学生运用龙贝格算法解决实际问题（塑料雨篷曲线满足函数y(x)=l sin(tx),则给定雨篷的长度后，求所需要平板材料的长度）。

二、方法原理2.1龙贝格积分原理龙贝格算法是由递推算法得来的。

由梯形公式得出辛普生公式得出柯特斯公式最后得到龙贝格公式。

在变步长的过程中探讨梯形法的计算规律。

设将求积区间[a ，b]分为n 个等分，则一共有n+1个等分点,k x a kh =+,0,1,b ah k n-==，n 。

这里用n T 表示复化梯形法求得的积分值，其下标n 表示等分数。

先考察下一个字段[1,k k x x +]，其中点()11212k k k xx x ++=+，在该子段上二分前后两个积分值()()112k k hT f x f x +=+⎡⎤⎣⎦ ()()21124k k k h T f x f x f x ++⎡⎤⎛⎫=++⎢⎥ ⎪⎢⎥⎝⎭⎣⎦显然有下列关系2112122k h T T f x +⎛⎫=+⎪⎝⎭将这一关系式关于k 从0到n-1累加求和，即可导出下列递推公式12102122n n n k k h T T fx -+=⎛⎫=+ ⎪⎝⎭∑ 需要强调指出的是，上式中的b ah n-=代表二分前的步长，而 1212k xa k h +⎛⎫=++ ⎪⎝⎭梯形法的算法简单，但精度低，收敛速度缓慢，如何提高收敛速度以节省计算量，自然式人们极为关心的。

根据梯形法的误差公式，积分值n T 的截断误差大致与2h 成正比，因此步长减半后误差将减至四分之一，既有21114n n T T -≈- 将上式移项整理，知2211()3n n n T T T -≈-由此可见，只要二分前后两个积分值n T 和2n T 相当接近，就可以保证计算保证结果计算结果2n T 的误差很小，这种直接用计算结果来估计误差的方法称作误差的事后估计法。

按上式，积分值2n T 的误差大致等于21()3n n T T -，如果用这个误差值作为2n T 的一种补偿，可以期望，所得的()222141333n n n n n T T T T T T =+-=- 应当是更好的结果。

2第二题 求数值积分11sin 3sin x dx x x -+⎰精确到610-。

龙贝格的算法思想：Romberg 方法也称为逐次分半加速法。

它是在复化梯形公式的基础上，利用Richardson 外推法，构造出一种高精度的数值积分方法。

在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行。

这样，前一次分割得到的函数值在分半以后仍可被利用,且易于编程 。

2.2.1程序设计关键步骤（1）准备初值aa[0][0]=(b-a)*(F(a)+F(b))/2.0来实现就算初值。

（2）用 for ( n=i,j=1; n>0;n--,j++)aa[j][n]=(pf(2*j+1)*aa[j-1][n+1]-aa[j-1][n])/(pf(j*2+1)-1) 来实现(2)(k )00T T 到的计算。

（3）用fabs(aa[i][1]-aa[i][0])判断数值的精度，如果fabs(aa[i][1]-aa[i][0]<e,则返回gcc(a,b,aa,i+1)进行计算下一个要用到得值； 如果fabs(aa[i][1]-aa[i][0]> e ，则根据aa[i+1][0]=(pf(2*i+3)*aa[i][1]-aa[i][0])/(pf(2*i+3)-1)而得到结果。

2.2.2程序运行结果如下图所示#include "stdafx.h"#include <iostream>#include <math.h>#define F(x) (sin(x)/(3*x+sin(x))) //函数举例。

using namespace std;//------------步长及4，16，64.......的实现----------double pf (int i){int s=1;for (int j=0;j<i;j++)s*=2;return s/2;}//---------定义一个求t1,t2......的函数-------------double gcc (double a, double b,double aa[][20],int i){double s,h,x;h=(b-a)/pf(i);s=0;x=a+h/3;do{s+=F(x);x+=h;}while (x<b);aa[0][i]=aa[0][i-1]/2+h/2*s;return 0;}//----------------------主函数---------------------------int main(){double aa[20][20]={0},e,a,b,h;int j,i,n;cout <<"请输入积分区间：\na= ";cin >>a;cout <<"b= ";cin >>b;cout <<"请输入精度：e=";cin >>e;aa[0][0]=(b-a)*(F(a)+F(b))/2.0;gcc(a,b,aa,1);aa[1][0]=(4*aa[0][1]-aa[0][0])/3;for (i=1;i<20;i++){gcc(a,b,aa,i+1);//求下一个要用的t。