数值积分_数值积分原理__matlab实现

课程设计报告课程设计题目：求解 的近似值课程名称：数值分析课程设计指导教师： X X X小组成员： X X XX X XX X X2013年12月31日目录目录 (1)题目 (2)一、摘要 (2)二、设计目的 (2)三、理论基础 (3)1、复合矩形法求定积分的原理 (3)2、复合梯形法求定积分的原理 (3)3、复合辛普森法求定积分的原理 (4)4、龙贝格求积公式原理 (5)四、程序代码及运算结果 (5)1、复合矩形法求定积分：用sum函数 (5)2、复合梯形法求定积分 (6)方法一 (6)方法二：用trapz函数 (7)3、复合辛普森法求定积分 (7)方法一 (7)方法二：用quad函数 (7)4、龙贝格求定积分 (8)5、Lobatto数值积分法 (9)6、波尔文(Borwein)高阶公式 (9)五、结果分析 (10)六、设计心得 (10)七、参考文献 (11)题 目：（1）已知：41102π=+⎰x dx ，所以 ⎰+=10214dx x π 。

于是，我们可以通过计算上述定积分的近似值来得到π的近似值。

（2）波尔文(Borwein )高阶公式在π值的高阶算法研究中，最好的结果来自两个都叫波尔文的数学家。

他们在1984年发表了一个2阶收敛公式：20=a ，00=b ，220+=p ，⎪⎪⎪⎩⎪⎪⎪⎨⎧++=++=+=++++++1111111)1()1(2)1(k k k k k kk k k k k k k b a b p p b a b a b a a a式中π→k p 。

试运用上述迭代算法，计算圆周率的近似值，并和前面传统方法进行比较。

一、摘要借助matlab 环境下的计算机编程语言，先用基本的积分函数对给出的题目进行求积分，然后基于给出的波尔文高阶收敛公式，在进行了连续迭代后，对运行结果做出分析，同时与之前的传统做法进行比较。

二、设计目的用熟悉的计算机语言编程，上机完成用复合矩形法、复合梯形法、复合辛普森法、龙贝格法以及Lobatto 数值积分方法，掌握各种方法的理论依据及求解思路，了解数值积分各种方法的异同与优缺点。

三、理论基础1、复合矩形法求定积分的原理定积分的几何意义是计算曲边梯形的面积，若将积分区间[a,b]分为n 段，每段都是一个小曲边梯形，用一个个小矩形代替这些小曲边梯形，然后把所有小矩形的面积加起来就近似得等于整个曲边梯形的面积，于是便求出了定积分的近似值，这就是复合矩形法求定积分的基本原理。

设第i 个小矩形的宽度),,2,1(1n i x x h i i i =-=-，高度)(i i x f y =，根据)()()()()()(00)0(00x f a b x f c a b dx x P dx x f baba-=-==⎰⎰可知其面积为i i i i i y h x f x x =--)()(1，可得出复合矩形求积公式：i ni i n n bah y h y h y h y f I dx x f ∑⎰==+++≈=12211)()(如果积分区间被等分为段：h n a b x x b x x x a i i n =-=-=<<<=-/)()(,110 ，则i ni i ni i bah y y n a b f I dx x f ∑∑⎰===-≈=11)()(。

2、复合梯形法求定积分的原理把积分区间[a,b]分成若干小区间，在每个小区间上以梯形面积近似曲边梯形面积，即用梯形公式求小区间上积分的近似值定积分存在定理表明，只要被积函数连续，当小区间长度趋于零时，小区间面积之和趋于曲边梯形面积的准确值，即定积分的准确值。

设第i 个小梯形的宽度),,2,1(1n i x x h i i i =-=-，两底高度分别为)(11--=i i x f y 和),,2,1)((n i x f y i i ==，则定积分的近似值为=+≈=∑∑⎰⎰=-=-i ni i i ni x x bah y y x f dx x f i i)(21)()(1111])()()()[(211112221110n n n n n n h y y h y y h y y h y y ++++++++---- 设被积函数)(x f 在],[b a 上连续可导，把],[b a 区间n 等分，令n a b h /)(-=，于是有⎪⎪⎭⎫ ⎝⎛++=+++++≈∑⎰-=-110121022)222(2)(n j j n n n bay y y h h y y y y y h dx x f 由于式中)(0a f y =，)(b f y n =，)()(jh a f x f y j j +==，代入上式得出复合梯形求积公式：n n j baT jh a f b f a f h dx x f =⎪⎪⎭⎫ ⎝⎛+++=∑⎰-=11)(2)()(2)( n T 表示区间分为n 等分时，用复合梯形求积法求出的定积分值。

3、复合辛普森法求定积分的原理 辛普森公式：)]()2(4)([6)(b f ba f a f ab f I +++-=因为辛普森公式用到了区间的中点，所以在构造复合辛普森公式时，把积分区间],[b a 等分为偶数份。

令m n 2=，其中m 为正整数，节点为)12,,2,1,0(-=+=m k kh a x k ，mab h 2-=在每两个小区间],[122+k k x x 上用辛普森公式，则有==∑⎰⎰-=+n m k xx bak kx f dx x f 102122)()(=+++∑∑-=+-=+12210122][)]()(4([3m k k k m k k k f R x f x f x f h∑∑∑-=-=-=++⎥⎦⎤⎢⎣⎡+++11010122][)(4)(2)()(3m k k m k m k k k f R x f x f b f a f h记⎥⎦⎤⎢⎣⎡+++=∑∑-=-=+10101222)(4)(2)()(3m k m k k k mx f x f b f a f h S 上式叫做复合辛普森公式，m S 2的下标m 2表示将积分区间],[b a m 2等分。

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

龙贝格秋季攻势是一个很实用的公式。

若已知n T 2与n T 的关系：)(21n n n H T T +=，其中∑-=-⋅=1112)(n i i n x f h H ，若记∑∑-=-=1112)(n i i xf 为新增分点函数值的和，则⋅=h H n ξ。

对n T 而言，∑为下标为奇数位置的函数值的和。

龙贝格求积公式：)()(0h T h T =，),2,1(12)(22)(2112 =--⎪⎭⎫⎝⎛=--k h T h T h T k k k k k 用⎪⎭⎫⎝⎛=i k i k h T T 2)(作为][f I 近似值，截断误差为)()1(2+k h O 。

四、程序代码及运算结果1、复合矩形法求定积分：用sum 函数编辑如下命令做出函数214x+，并保存： function y= qiupai(x); y=4./(1+x.^2);该函数的图像如右图所示，生成方法如下：x=0:0.01:1;plot(x,qiupai1(x),'linewidth',2), gridlegend('qiupai1(x)')在命令窗口输入并运行：h=1/100; x=0+h:h:1;s=h*sum(qiupai1(x));pi1=vpa(s,6); format longpi1运行结果：pi1 = 3.13158当将分割的小矩形的宽度变小时, π值的精度相应提高,下表即为改变h的值得到的相应的pil的值：2、复合梯形法求定积分方法一先编辑梯形函数并保存：function s=tixing(f,a,b,n)% f是被积函数；% a,b分别为积分的上下限；% n是子区间的个数；% s是梯形总面积；h=(b-a)/n;s=0;for k=1:(n-1)x=a+h*k;s=s+feval('f',x);endformat longs=h*(feval('f',a)+feval('f',b))/2+h*s;然后在编辑窗口输入如下命令并运行：pi1=tixing('4/(1+x^2)',0,1,100)运行结果：pi1 = 3.141575986923129当将分割的小矩形的宽度变小时, π值的精度相应提高,下表即为改变h的值得到的相应的pil的值：方法二：用trapz函数在命令窗口输入：x=0:0.01:1;format longpi2=trapz(x,qiupai1(x))% 函数qiupai1(x)已在前面给出并保存运行结果：pi2 = 3.1415759869231283、复合辛普森法求定积分方法一先编辑辛普森函数并保存：function s=simpson(f,a,b,n)% f是被积函数；% a,b分别为积分的上下限；% n是子区间的个数；% s是梯形总面积,即所求积分数值；h=(b-a)/(2*n);s1=0;s2=0;for k=1:nx=a+h*(2*k-1);s1=s1+feval('f',x);endfor k=1:(n-1)x=a+h*2*k;s2=s2+feval('f',x);ends=h*(feval('f',a)+feval('f',b)+4*s1+2*s2)/3;然后在编辑窗口输入如下命令并运行：pi3 = simpson('4/(1+x^2)',0,1,100)运行结果：pil3 = 3.141592653589793方法二：用quad函数求qiupai1(x)从0~1积分，可编辑如下的命令：format long;pi3=quad('qiupai1',0,1)运行结果：pi3 = 3.1415926829245674、龙贝格求定积分先编辑如下命令并保存：function [R,pil4,err,h]=romber(f,a,b,n,delta)% f是被积函数% a,b分别是积分的上下限% n+1是T数表的列数% delta是允许误差% R是T数表% pil4是所求积分值M=1;h=b-a;err=1J=0;R=zeros(4,4);R(1,1)=h*(feval('f',a)+feval('f',b))/2while ((err>delta)&(J<n))|(J<4)J=J+1;h=h/2;s=0;for p=1:Mx=a+h*(2*p-1);s=s+feval('f',x);endR(J+1,1)=R(J,1)/2+h*s;M=2*M;for K=1:JR(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K))/(4^K-1);enderr=abs(R(J,J)-R(J+1,K+1));endpil4=R(J+1,J+1)然后在编辑窗口输入如下命令并运行：romber('4/(1+x^2)',0,1,100,0.5*(10^(-8)))运行结果：pil4 = 3.141592653589723ans的值是如下的矩阵：3 0 0 0 0 0 03.1 3.133333 0 0 0 0 0 3.131176 3.141569 3.142118 0 0 0 0 3.138988 3.141593 3.141594 3.141586 0 0 0 3.140942 3.141593 3.141593 3.141593 3.141593 0 0 3.14143 3.141593 3.141593 3.141593 3.141593 3.141593 0 3.141552 3.141593 3.141593 3.141593 3.141593 3.141593 3.1415935、Lobatto数值积分法在要求的绝对误差范围内，用自适应递推步长复合Lobatto数值积分法，与它相应的是高阶数值积分函数。