数值分析实验 复化求积公式(改进版)
复化求积公式
h[ 2
f ( x0 )
n1
2
k 1
f ( xk )
f ( xn )]
复化梯形公式
Tn
h 2
[
f
(
x0
)
n1
2
k 1
f ( xk )
f ( xn )]
复化梯形公式
计算方法
2.复化辛浦生公式
计算方法
在每个小区间[xk1, xk ]上应用辛浦生公式得：
xk
xk 1
则
f
b
( x)dx h[ f 6
计算方法
在 每 个 小 区 间[ xk1, xk ]上 应 用 梯 形 公 式 得 ：
xk xk 1
f ( x)dx
h 2
[
f
(
xk1
)
f ( xk )]
则
b
n
f (x)dx =
xk f (x)dx
a
k 1 xk1
n k 1
h[ f 2
(xk1)
f
(xk )]
计算方法
x0 x1 x2 x3
2
三、区间逐次分半求积法
计算方法
复化求积公式可有效提高计算精度，但对给定 的误差限，如何确定节点的个数，即[a,b]应多少等 份？由截断误差可以估计步长的取值情况，但需要 给出各阶导数的最大值，这往往是比较困难的，且 估计值往往偏大.
接下来，我们将考虑步长的更为实用的选取方 法.
计算方法
若用Tn及T2n分别表示将[a, b]n等分及2n等分的复化 梯形公式，则
f(x) 1 0.997 0.9896 0.976 0.95 0.936 0.908 0.877 0.841 3978 158 7267 8851 1556 8516 1925 4709
chap4第2节 复化求积公式
Rn [ f ]
h (b a )
2
f ( ), (a , b)
12
如果记 M 2 max f ( x )
a xb
则有 Rn [ f ]
b
a
f ( x )dx Tn
( b a )h 12
2
M2
(b a ) 12n
2
3
M2
上式说明复化梯形公式是收敛的。
这时由
xk x k 1
得到
h h f f ( x )dx f ( xk 1 ) 4 f ( x 1 ) f ( xk ) k 2880 6 2 n
5
(4)
( k )
b
a
f ( x )dx
i 1
xk
f ( x )dx
x k 1
5 n h h (4) f ( k ) f ( x k 1 ) 4 f ( x 1 ) f ( x k ) 2880 k 1 k k 1 6 2
1
1 4 4 4 2 2 2 1 4 6 1 1 9 9
3.1230
4 )3
而梯形公式的结果为
1 x
0
1
4
2
dx
1 0 2
(
4
1 0 11
例 4.4 用复化梯形公式计算积分 I 0 e dx ，应将区间 [0,1]多少等分，才可以使其截断误差不超过 1 10 4
x
1
2
解：复化梯形公式的误差为
Rn [ f ] f ( x )dx Tn
a b
(b a ) 12n
复化求积公式的算法及其应用
摘要在数值计算中，低阶牛顿柯特斯求积方法存在很多缺陷，从余项公式可以看出其要求提高求积公式的代数精度，必须增加结点个数，会导致插值多项式出现龙格现象，且数值稳定性不能保证.基于以上原因，我们往往采用复化求积方法，此方法不仅可以克服以上缺点而且便于在计算机上实现，值得研究和学习.在本课程设计中，我们首先从复化求积公式的思想引入，然后详细介绍复化梯形求积公式、复化辛普森求积公式和复化柯特斯求积公式的推导过程和相关性质，再对三种求积公式进行比较和总结，其次画出三种求积公式的流程图，最后通过求解例题写出三种求积算法的程序设计.关键词复化求积算法;流程图;程序设计目录引言 0第一章复化求积算法 (2)§1.1复化求积公式 (2)§1.1复化求积公式的思想 (3)§1.2复化求积公式的构造 (3)§1.2复化梯形求积公式 (3)§1.2.1复化梯形求积公式的推导过程 (3)§1.2.2复化梯形求积公式的性质 (3)§1.3复化辛普森求积公式 (4)§1.3.1复化辛普森求积公式的推导过程 (4)§1.3.2复化辛普森求积公式的性质 (4)§1.4复化柯特斯求积公式 (5)§1.4.1复化柯特斯求积公式的推导过程 (5)§1.4.2复化柯特斯求积公式的性质 (5)§1.5三种复化求积公式的比较及总结 (6)第二章复化求积公式算法的流程图及其应用 (9)§2.1 流程图 (9)§2.2 应用 (12)参考文献 (15)附录A (16)附录B (17)附录C (18)引言积分计算在分析数学领域里是个古老的问题，在数值分析中已被广泛应用.但在计算机上却不能像在分析数学中那样，用原函数[满足)()('x f x F =的函数)(x F 就是函数)(x f 的原函数]计算积分.这是因为在实际问题中，函数关系往往是用列表数据或曲线给出的.即使知道了函数的表达式，求其一个原函数并非一个简单问题.许多函数难以用初等函数表示（如2,/sin x e x x -等）.在计算机上，通常利用函数的若干个离散值，以代数运算近似计算积分值，这类近似计算法称为数值积分法.设给定区间],[b a 上的函数)(x f .需要建立计算积分dx x f f I ba ⎰=)()(的近似方法.数值积分的基本思想是试图用一个简单又易于积分的函数逼近)(x f ，以计算积分)(f I .显然插值多项式是一个很好的选择，因为插值多项式可由)(x f 的若干值构造出来，其积分很容易计算.为此，需将],[b a 分为n 等分n i x x i i ,,2,1],,[1 =+，其中b x x x x a n =<<<<=+1321 .分割步长h ，因此,1,3,2,/)1(1+=-+=n i h i x x i 对应的函数值)()(,),(),()(121b f x f x f x f a f n ==+ .显然)(f I 可以表示为所有小区间上各函数的积分的和，即)()(1f I f I ni i ∑==其中 dx x f I i ix x i ⎰+=1)(通常把为每个)(f I i 建立的计算公式简称为求积公式，而把)(f I 建立的求积公式称为复化求积公式.由于在实际计算时，不宜使用高阶的牛顿——柯特斯公式，但若积分区间较大，单独用一个低阶的牛顿——柯特斯公式来计算积分的近似值，显然精度不好，为了提高数值求积的精确度，可利用积分对区间的可加性来解决这个问题，这就是通常采用的复合求积法.而且使用这种方法之后，求积公式的收敛性和稳定性也得到了改善.第一章 复化求积算法牛顿—柯特斯公式的求积余项表明，求积节点n 越大，对应的求积公式精度越高，但由于牛顿—柯特斯公式在8>n 时数值不稳定，因此不能用增加求积节点数的方法来提高计算精度.实用中常将求积区间],[b a 分成若干个小区间，然后在每个小区间上采用数值稳定的牛顿—柯特斯公式求小区间上的定积分，最后把所有小区间上的计算结果相加来作为原定积分的近似值.采用这种方法构造的求积公式就称为复合求积公式.复合求积公式具有计算简单且可以任意逼近所求定积分值的特点，这是牛顿—柯特斯公式一般做不到的.常用的复合求积公式有复合梯形求积公式和复合辛普森求积公式以及复合柯特斯求积公式.以下我们将从三种复化求积算法的构造、余项、稳定性、收敛性等几方面进行讨论，并写出相应的流程图以及应用中所涉及到的算法的程序设计.§1.1复化求积公式§1.1.1 复化求积公式的思想n 很大时，牛顿——柯特斯求积公式出现了不稳定、不收敛现象，往往使用低阶牛顿——柯特斯求积公式，误差比较大，故将],[b a 若干等分，在每个子区间上反复使用低阶牛顿——柯特斯公式，进行累加.而构造出来的新的求积公式，称之为复化求积公式.在构造求积公式的过程中，我们将求积区间],[b a 进行等距细分：n i nab ia x i ,,1,0, =-+=，在每个小区间],[1i i x x -上用相同的“基本”求积公式(如梯形公式；中矩形公式；左(右)矩形公式或辛普森公式)计算出dx x f i i x x ⎰-1)(的近似值i S .§1.1.2 复化求积公式的的构造将定积分⎰ba dx x f )(的区间],[b a 划分为n 等分，各节点为kh a x k +=，n k ,,1,0 =，nab h -=，在子区间)1,,1,0](,[1-=+n k x x k k 上使用牛顿——柯特公式，将],[1+k k x x 分割为l 等份，步长为l h，节点为1,,2,,+=+++k k k k k x llhx l h x l h x x记121,,,,++++=k ll k lk lk k x xxxx为在],[1+k k x x 上作)(x f 的l 阶牛顿——柯特斯求积公式.∑∑⎰=++=+=-=≈+li li k l i li k li l i k k k i x x xf C h xf C x x I dx x f k)(0)(1)()()()()(1由积分区间的可加性，可得nli k n k li l i n k k l n k k k baI xf C h I dxx f dx x f ==≈=+-==-=-=+∑∑∑∑⎰⎰)()()(100)(1)(11§1.2 复化梯形求积公式§1.2.1 复化梯形求积公式的的推导过程将积分区间],[b a 划分等分，步长nab h -=，求积节点kh a x k +=，n k ,,1,0 =在每个小区间)1,,1,0](,[1-=+n k x x k k 上应用梯形公式)]()([2)(11++≈⎰+k k x x x f x f hdx x f k k然后将它们累加求和，作为所求积分I 的近似值.])()(2)([2)]())()()((2)([2)]()([2)()(11121011011∑∑∑⎰⎰---+-=-=++=+++++=+≈==+n i k n n k k n k n k x x bab f x f a f hx f x f x f x f x f hx f x f hdx x f dx x f I k k记n T )]()(2)([211b f x f a f hn i k ++=∑-=式为复化梯形求积公式，下标n 表示将区间n 等分，若把区间n 2等分，在每个小区间上仍用梯形求积公式，则可得到n n T T ,2和n H 间的关系为：)(212n n n H T T +=其中∑=--+=nk n nab k a f h H 1]2)12([ §1.2.2复化梯形求积公式的性质性质1.1复化梯形求积公式余项当)(x f 在],[b a 上有连续的二阶导数，则复化梯形公式的余项：)(12)()(''2ηf h a b T dx x f R n ba T --=-=⎰ ],[b a ∈η 性质1.2稳定性若],[,)(''b a x M x f ∈≤，则有估计式M na b R nT 2312)(-≤ 复化梯形求积公式的系数均大于零，且满足a b nh n hA ni i -==+-+=∑=]1)1(21[2因此，复化梯形求积公式的计算过程是数值稳定的.性质1.3收敛性可证复化梯形求积公式是收敛的. 性质1.4代数精度定义1.1 若积分⎰b adx x f )(的数值积分公式⎰badx x f )()(0k nk k x f A ∑=≈对于任意一个次数不高于m 次的多项式都精确成立，且存在一个1+m 次多项式使之不精确成立，则称该数值积分公式的代数精度为m .可证复化梯形求积公式的代数精度为2.§1.3 复化辛普森求积公式§1.3.1 复化辛普森求积公式的的推导过程将积分区间],[b a 划分等分，记子区间],[1+k k x x 的中点为h x x k k 2121+=+在每个小区间上应用辛普森公式，则有))()(2)(4)((6)444(6)]()(4)([6)()(101121211223112101211011b f x f x f a f hf f f f f f f f f hx f x f x f hdxx f dx x f I n k n k k k n n n k k k n k n k x x bak k+++=+++++++++=++≈==∑∑∑∑⎰⎰-=-=+--++-=-=+其中h x xk k 2121+=+记 )]()(2)(4)([6111021b f x f x f a f hS n k k n k k n +++=∑∑-=-=+式为复化辛普森求积公式§1.3.2复化辛普森求积公式的性质性质1.5复化辛普森求积公式余项当)(x f 在],[b a 上有连续的四阶导数，复化辛普森公式的求积余项为：)(2880)()2(180)4(4)4(4ηηf h a b f h a b R S --=--= ],[b a ∈η 性质1.6稳定性同复化梯形求积公式，复化辛普森求积公式的系数均大于零，且满足总和为a b - 因此，复化辛普森求积公式的计算过程是数值稳定的.性质1.7收敛性可证复化辛普森求积公式是收敛的. 性质1.8代数精度可证复化辛普森求积公式的代数精度为4.§1.4 复化柯特斯求积公式§1.4.1 复化柯特斯求积公式的的推导过程将积分区间],[b a 划分等分，若把每个子区间],[1+k k x x 四等份，内点依次记为432141,,+++k k k xxx，同理可得复化柯特斯求积公式)](7)(14)(32)(12)(32)(7[9010101143211041b f x f x f x f x f a f hC n k n k n k k k k n k k n +++++=∑∑∑∑-=-=-=++-=+(1-1)其中h x xh x x h x x k k k k k k 43;21;41432141+=+=+=+++ 记(1-1)为复化柯特斯求积公式§1.4.2复化柯特斯求积公式的性质性质1.9复化柯特斯求积公式余项当)(x f 在],[b a 上有连续的四阶导数，复化柯特斯公式的求积余项为：)()4(945)(2)6(6ηf h a b R c --= ],[b a ∈η性质1.10稳定性同复化梯形求积公式，复化柯特斯求积公式的系数均大于零，且满足总和为a b - 因此，复化柯特斯求积公式的计算过程是数值稳定的.性质1.11收敛性可证复化柯特斯求积公式是收敛的. 性质1.12代数精度可证复化柯特斯求积公式的代数精度为6.§1.5 三种复化求积公式的比较及总结为了更形象的表述三种复化求积公式之间的关系，我们通过一个例子来进行比较例1.1使用各种复化求积公式计算定积分dx xxI ⎰=10sin 为简单起见，依次使用8阶复化梯形公式、4阶复化辛普森公式和2阶复化柯特斯公式,可得各节点的值如下表表1-1节点值94569086.0)]1()(2)0([161718=++=∑=f x f f T k k 94608331.0)]1()(2)(4)0([2413031214=+++=∑∑==+f x f x f f S k k k k 94608307.0)]1(7)(14)](32)(12)(32[)0(7[180111104342412=+++++=∑∑==+++f x f x f x f x f f C k k k k k k 比较三个公式的结果：精度最低 94569086.08=T 精度次高 94608331.04=S 精度最高 94608307.02=C原积分的精确值为6719460830703.0sin 10==⎰dx xxI . 我们知道，三种求积公式的余项分别如表1-2表1-2 复化梯形、辛普森、柯特斯求积公式的余项定义1.2对于复化求积公式n I 若存在0>p 及0≠c ，使其余项n I I -满足c h I I pnh =-→0lim则称复化求积公式n I 是p 阶收敛的 P 阶收敛性的意义：对于一个数值求积公式来说，收敛阶越高，近似值n I 收敛到真值dx x f ba ⎰)(的速度就越快.由于三种求积公式的余项分别是h 的2,4,6阶无穷小量 所以n n n C S T ,,趋于定积分I 的速度依次更快.从这三种求积公式的构造过程中可以看出，它们都属于机械求积公式，但不属于插值行和牛顿柯特斯公式.都具有稳定性和收敛性，且收敛速度一个比一个快，一个比一准确.在使用函数值个数相等的情况下，248,,C S T 的精度逐渐升高.第二章复化求积公式算法的流程图及其应用§2.1 流程图1.复化梯形求积公式图2.1 复化梯形求积公式算法的流程图Step1给出被积函数)(x f 、区间],[b a 端点b a ,和等分数n ； Step2求出,kh x k =nab h -=； Step3计算∑-=1)(),(),(n k k x f b f a f ；Step4得)]()()([211b f x f a f h T n k k n ++=∑-=2. 复化辛普森求积公式图2.2 复化辛普森求积公式算法的流程图Step1 给出被积函数)(x f 、区间],[b a 端点b a ,和等分数n ； Step2求出,kh x k =nab h -=； Step3计算∑∑-=+-=1211)(,)(),(),(n k k n k k xf x f b f a f ；Step4得)]()(2)(4)([6111021b f x f x f a f hS n k k n k k n +++=∑∑-=-=+3. 复化柯特斯求积公式图2.3 复化柯特斯求积公式算法的流程图Step1给出被积函数)(x f 、区间],[b a 端点b a ,和等分数n ； Step2求出,kh x k =nab h -=； Step3计算∑∑∑∑-=-=+-=+-=+11143121141)(,)(,)(,)(),(),(n k k n k k n k k n k k x f xf xf xf b f a f ；Step4得)](7)(14)(32)(12)(32)(7[9010101143211041b f x f x f x f x f a f hC n k n k n k k k k n k k n +++++=∑∑∑∑-=-=-=++-=+§2.2 应用例2.1.分别用复化梯形，复化辛普森，复化柯特斯公式计算函数32)(x x x f -=在区间]1,0[上的弧长S .(要求写出源程序和运行结果) *注 在],[b a 上的弧长dx x f S ba⎰+=2'))((11.用复化梯形公式计算S 的过程：(1).写出变量说明表2-1 复化梯形求积公式程序设计的变量说明Step1 输入n ，nab h -=，被积函数0),(1=s x f ； Step2 for 1=k to 1-n ；{计算11)(s kh a f s →++}))(2)((21b f s a f hs ++=； Step3 输出近似值s .(3) 写出源程序和运行结果(见附录A) 2.用复化辛普森公式计算S 的过程: (1).写出变量说明表2-2 复化辛普森求积公式程序设计的变量说明Step1:输入n ，nab h -=，被积函数0),(1=s x f 0,2=s ； Step2:for 1=i to 1-n ，2+=i i ；{计算11)2/*(s h i a f s →++} Step3:for 2=j to 1-n ，2+=j j ； {计算22)2/*(s h j a f s →++}))(24)((621b f s s a f hs +++=； Step4:输出近似值s .(3).写出源程序和运行结果(见附录B) 3.用复化柯特斯公式计算S 的过程: (1).写出变量说明表2-3 复化柯特斯求积公式程序设计的变量说明Step1输入n ，nab h -=，被积函数0),(1=s x f 0,2=s 0,3=s ； Step2 for 1=i to 1-n ，2+=i i ；{计算11)4/*(s h i a f s →++} Step3:for 2=j to 1-n ，4+=j j ； {计算22)4/*(s h j a f s →++} Step4: for 4=k to 2-n ，2+=k k ； {计算33)4/*(s h k a f s →++}))(141232)((90321b f s s s a f hs ++++=； Step5:输出近似值s .(3).写出源程序和运行结果(见附录C)根据运行结果可知，由三种复化求积公式求得的S 的值分别为064837.1、061199.1、061189.1，精度逐渐升高.参考文献[1] 薛毅，耿美英.数值分析[M]. 北京:北京工业大学出版社.2003年. [2] 刘长安.数值分析教程[M].西安:西北工业大学出版社.2005年.[3] 朝伦巴根，贾德彬.数值计算方法[M].北京:中国水利水电出版社.2007年.[4] 韩旭里，万中.数值分析与实验[M].北京: 科学出版社.2006年.[5] 林成森.数值分析[M].北京: 科学出版社.2007年.[6] 封建湖，车刚明，聂玉峰.数值分析原理. 北京: 科学出版社.2001年.附录A1.复化梯形求积公式的程序设计：(1).源程序：#include<stdio.h>#include<math.h>double f(double x){double z;z=sqrt(1+pow((2*x-3*pow(x,2)),2));return z;}main(){ int n,k;float h;float a;float b;double s=0.0;double s1=0.0;double t;printf("Please input the deng fen ;");scanf("%d",&n);printf("Please input qujian a ;");scanf("%f",&a);printf("Please input qujian b ;");scanf("%f",&b);h=(b-a)/n;for (k=1;k<n;k++){ t=a+k*h;s1=s1+f(t);}s=(h/2)*(f(a)+2*s1+f(b));printf("%f\n",s);}(2).运行结果：图1 复化梯形求积公式计算弧长结果附录B2.复化辛普森求积公式的程序设计：(1).源程序：#include<stdio.h>#include<math.h>double f(double x){double z;z=sqrt(1+pow((2*x-3*pow(x,2)),2));return z;}main(){ int n,i,j;float h;float a;float b;double s=0.0;double s1=0.0,s2=0.0;double t,l;printf("Please input the deng fen ;");scanf("%d",&n);printf("Please input qujian a ;");scanf("%f",&a);printf("Please input qujian b ;");scanf("%f",&b);h=(b-a)/n;for(i=1;i<8;i=i+2){t=a+i*h/2;s1=s1+4*f(t);}for(j=2;j<8;j=j+2){l=a+j*h/2;s2=s2+2*f(l);}s=(h/6)*(f(a)+s1+s2+f(b));printf("%f\n",s);}(2).运行结果：图2 复化辛普森求积公式计算弧长结果附录C3.复化柯特斯求积公式的程序设计：(1).源程序：#include<stdio.h>#include<math.h>double f(double x){double z;z=sqrt(1+pow((2*x-3*pow(x,2)),2));return z;}main(){int n,i,j,k;float h;float a;float b;double s=0.0;double s1=0.0,s2=0.0,s3=0.0;double t,l,m;printf("Please input the deng fen ;"); scanf("%d",&n);printf("Please input qujian a ;");scanf("%f",&a);printf("Please input qujian b ;");scanf("%f",&b);h=(b-a)/n;for(i=1;i<8;i=i+2){t=a+i*h/4;s1=s1+32*f(t);}for(j=2;j<7;j=j+4){l=a+j*h/4;s2=s2+12*f(l);}for(k=4;k<6;k=k+2){m=a+k*h/4;s3=s3+14*f(m);}s=(h/90)*(7*f(a)+s1+s2+s3+7*f(b)); printf("%f\n",s);}(2).运行结果：图3 复化柯特斯求积公式计算弧长结果。
数值分析63 复化求积公式龙贝格求积公式讲解
只增加了一个分点
1 xk?1/ 2 ? 2 ( xk ? xk?1)
设hn=(b? a)/n, xk=a+kh n (k=0,1,? ,n),在[xk, xk+1] 上用梯形公式得
T1 ?
hn 2
?f
(
xk
)
?
f ? ( xk ? 1 )
复化求积的基本想法 :
将积分区间 [a, b]n等分, 步长
h?
b
? n
a
,
分点为
xk=a+kh (k=0,1,…,n) , 则由定积分性质知
? ? ? I ?
b
n?1
f ( x )dx ?
xk?1 f ( x )d x
a
k ? 0 xk
每个子区间 上的积分
?xk?1 f ( x )dx xk
用低阶求积公式 , 然后把所有区间的 计算结果求和 ,
注2: 同样也可用 | S4m-S2m |<ε 来控制计算的精度 . 这就是下面要介绍的 龙贝格求 积公式 .
6.4 龙贝格求积公式
6.4.1 梯形公式的递推化
复化求积方法可提高求积精度，实际计算时若
精度不够可将步长逐次分半 . 设将区间 [a, b]分为n等
分，共有 n+1个分点，如果将求积区间再分一次，则 分点增至 2n+1个，我们将二分 前后两个积分值 联系
果T8=0.9456909 只有2位有效数字，而应用复化辛普 森公式计算的结果 S4= 0.9460832 却有6位有效数字 .
注：为了利用余项公式估计误差，要求 f(x)=sin x/x 的高阶导数，由于
数值分析63 复化求积公式龙贝格求积公式讲解
将积分区间 [a, b] 划分为2n等分, 即将每一个区间
[xk, xk+1]经过二等分增加了一个分点
x k ?1/2
?
1 2
(
x
k
?
x k?1)
在每个子区间 [xk, xk+1]上的积分用 辛普森公式 , 得
?x k ? 1 xk
f (x)dx ?
h 6
[
f
(
xk
)
?
4
f
(
xk
? 1/2
)
?
f (xk?1)]
?
(b ? a)5 2880n4
f
(4) (? )
?6.3.2 复化辛普森公式
将积分区间 [a, b] 划分为2n等分, 则
? ? ? I ?
b
n?1
f ( x )dx ?
x2k? 2 f ( x )d x
a
k ? 0 x2k
每个子区间 [x2k, x2k+2]上的积分用 辛普森公式 , 得
?x2 k ? 2 x2k
称为复化辛普森公式 . 记
? ? h
n?1
n?1
Sn ?
[f 6
(a) ?
4
k?0
f
( x k ? 1/2 ) ?
2
k ?1
f
(xk ) ?
f (b)]
若 f(x)? C 4[a,b], 其求积余项 为 h ? b ? a
n
Rn ( f ?
b ? a ( h )4 180 2
f (4) (? ) ?
f (b)]
称为复化梯形公式 . 记
? h
n?1
Tn
?
复化梯形公式和复化辛普生公式
void main()
{
tixing ti;
ti.integration();
}
2）复化simpson
#include<iostream.h>
#include<math.h>
class simpson
{
private:
double result;
int n,k;//n是等分数，
double h;//h步长
cin>>b>>a;
cout<<"输入你要使用tixing法则的数目(即等分数)";
cin>>n;
h=(b-a)/n;
sum=0.0;
for(k=1;k<n;k++)
{
sum+=func(k*h);
}
integral=(h/2.0)*(2.0*sum+func(b)+1);
cout<<endl<<"积分值="<<integral<<endl;
public:
double sine(double);
void integration();//实现积分的函数
};
double simpson::sine(double a)
{
result=sin(a)/a;//被积函数
cout<<"sin("<<a<<")/"<<a<<"="<<result<<endl;
源程序：
复化求积公式
2 定理7 定理 若 f ( x ) ∈ C [a , b ] , 则复化梯形公式的余项为
说明： ） （ 式说明复化梯形公式的余项收敛于零的速度与 说明： 1）(3.3)式说明复化梯形公式的余项收敛于零的速度与 h2收敛于零的速度相同，即余项等 于O(h2)。 收敛于零的速度相同， 。 导数值)确定 （2）余项可由端点的函数值 导数值 确定。 ）余项可由端点的函数值(导数值 确定。 3、复化中矩求积公式 （推导类似复化梯形公式） 、 推导类似复化梯形公式） x 1 x 1 x3 上采用中矩形公式， 在 [ x i −1 , x i ]上采用中矩形公式， n− 2 2 2 b − a xi xi −1 + xi , 记h = xn =b x = ∫xi−1 f ( x)dx ≈ h ⋅ f ( 2 )， a=x0 x1 x2 n i = 1,2,L, n 所以
即
∫
b
a
n −1 x i ) + f (b) 2 i =1 2
记
= Tn
( 3 .1)
下面考虑余项，先从每个小区间上考虑余项， 下面考虑余项，先从每个小区间上考虑余项，因为每个小区 间上是N－ 公式中当 公式中当n=1时的梯形公式。 时的梯形公式。 间上是 －C公式中当 时的梯形公式
n n
所以 − h∑ f ′′(ξ i ) = −∑ O( h2 ) + f ′(a ) − f ′(b) → f ′(a) − f ′(b), (h → 0)，
n
i =1 i =1
即 − h∑ f ′′(ξi ) → f ′(a) − f ′(b) (h →0)。
i=1
#
− (b − a)h2 ∫a f ( x)dxb − Tn = 12 f ′′(ξ ), a < ξ < b， ( 3 . 2 ) ∫a f ( x)dx − Tn → 1 ( f ′(a) − f ′(b)), h → 0。 及渐近估计式 ( 3 .3 ) 2 12 h
复化求积公式求数值积分
应用软件课程设计——复化求积公式求数值积分数学0801班 ***一、实验目的程序1：利用复化梯形公式以及复化辛普森公式求解定积分的数值解。
程序2：分析剖分区间个数对复化梯形公式精度的影响。
程序3：比较MATLAB 的quad 、quadl 命令与上述两种方法的精度；比较在相同剖分区间下两种求数值积分方法的精度。
分析与探讨两种方法精度不同的原因。
说明：原题目给出的积分精确值I=4.006994稍过粗糙，所以通过计算机求解得到更为精确的解析解。
详见测试结果。
二、算法说明自定义函数有：积分函数〔在此题中为dx e I x ⎰+=201〕：hanshu.m 〔附录1〕复化梯形公式求定积分数值解tixing.m 〔附录2〕 复化辛普森公式求定积分数值解xps.m 〔附录3〕 积分原函数解析解jxj.m 〔附录4〕程序1：p1.m 〔附录5〕1、n=2000;%确定积分区间分割份数2、tx=tixing(0,2,n);%用复化梯形公式求解3、xps=xps(0,2,n);%用复化辛普森公式求解4、显示结果disp(['积分区间分割分数为：',num2str(n)]) disp('复化梯形公式的求解结果：'),disp(tx) disp('复化辛普森公式的求解结果：'),disp(xps)程序2-1：p21.m 〔附录6〕1、jxj=jxj(2)-jxj(0);%求出解析解。
2、tx=zeros(5,1);%给数值解向量赋值。
d=zeros(5,1);%给误差向量赋值。
3、i=1,2,3,4,54、n=10^i;%定义剖分份数。
5、tx(i)=tixing(0,2,n);%将剖分份数n 代入，求出该n 下的数值解。
d(i)=tx(i)-jxj;%求出误差。
6、结束循环 7、输出结果程序2-2：p22.m 〔附录7〕1、jxj=jxj(2)-jxj(0);%求出解析解。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
-2 x2 1 f 3 3x f1
3
f2
4 1 x2 f 4 xe x
(1) ln 2 ln 3 2
因为
f1
(2)
3
2
1 dx x2 1
-4(3x 2 1) ( x 2 1)3
f 1(4) (
24 24 ) 5 ( x 1) ( x 1)5
所以
max f 1(2) = f 1(2)
x 2
52 27 1 2
max f 1(4) = f 1(4)
x 2
5808 243
又因为题目要求绝对误差限为 107 ，所以 对于复化梯形求积公式有
(b a) 2 '' (3 2)h2 52 1 Rn ( f ) h f ( ) = 10-7 12 12 27 2
1 2 (3) 3x dx ln 3 0
ba 1 n 26
因为
f 3(2) 3x (ln3)2
所以
max( f 3(2) ) f 3(2)
x 1
f 3(4) 3x (ln 3)4
3(ln 3)2
1 2
max( f 3(4) ) f 3(4)
x 1
3(ln 3)4
n 1
2.复化 Simpson 求积公式
将区间[a,b]划分为n等分，在每个子区间[xk , xk 1 ]上采用Simpson式 f ( x)dx
a b
(b a) ab 1 [ f (a) 4 f ( ) f (b)], 若记xk 1 2 =xk h,则得 6 2 2 n 1 n 1 b xk 1 h I f ( x) dx f ( x) dx [ f ( xk ) 4 f ( xk 1 2 ) f ( xk 1 )] Rn ( f ) a x 6 k 0 k 0 k 记
所以
n 18.2
因此 取节点数 n=19 步长 h
ba 1 n 19
1 (2) 4 dx 0 1 x2
1
因为
f 2(2) 8(3x 2 1) ( x 2 1)3 f 2(4) 96(5 x 4 10 x 2 1) ( x 1)5
所以
max f 2(2) = f 2(2)
所以
n 1791.6
因此 取节点数 n=1792 步长 h 对于复化 Simpson 求积公式有
Rn ( f ) b a h 4 (4) (3 2)h 4 5808 1 ( ) f ( ) = 10-7 4 180 2 180 2 243 2
ba 1 n 1792
所以
n 20.1
因此 取节点数 n=21 步长 h
ba 1 n 21
4
对于复化 Gauss-Legendre 型求积公式有
Rn ( f )
(b a)h 4 (4) (3 2)h 4 5808 1 f ( ) = 10-7 4320 4320 243 2
Tn =
b 1 1 (b a )3 '' R[f ]= f ( x) dx (b a)[ f ( a) f (b)] f ( ) (a b) a 2 2 12
得 h3 '' f (k ) ], k [xk , xk 1 ] 12 k 0 1 n 1 由于f ( x) C 2 [a, b], 且 min f '' (k ) f '' (k ) max f '' (k ) 0 k n 1 0 k n 1 n k 0 所以 (a, b), 使 Rn ( f )= I -Tn = [ f '' ( )= 于是复化梯形公式余项为 Rn ( f )= (b a ) 2 '' h f ( ) 12 1 n 1 '' f (k ) n k 0
ba 1 n 14
对于复化 Gauss-Legendre 型求积公式有
Rn ( f ) (b a)h 4 (4) (1 0)h 4 1 f ( ) 3(ln 3)4 = 10-7 4 4 2 270 270 2 2
所以
n 11.93
因此 取节点数 n=12 步长 h
ba 1 n 24
8
对于复化 Gauss-Legendre 型求积公式有
Rn ( f ) (b a)h 4 (4) (1 0)h 4 2 1 f ( ) 6e = 10-7 4 4 2 270 270 2 2
Rn ( f ) b a h 4 (4) (1 0)h 4 1 ( ) f ( ) 3(ln 3)4 = 10-7 4 180 2 180 2 2
ba 1 n 2457
所以
h ba 1 7.58 10-2 n n n 13.19
因此 取节点数 n=14 步长 h
二、复化求积公式求解过程
1.利用余项对所要求的每种算法做出步长的事前估计
(1) ln 2 ln 3 2 (3)
1 2 3x dx ln 3 0
3
2
1 1 1 d x (2) 4 dx 2 2 0 x 1 1 x
(4) e2 xe x dx
1
2
根据题意： 令
ba 1 n 12
(4) e xe x dx
2 1
7
2
因为
f 4(2) 2e x xe x f 4(4) 4e x xe x
所以
max( f 4(2) ) f 4(2)
x2
4e2
1 2
max( f 4(4) ) f 4(4)
x2
6e2
又因为题目要求绝对误差限为 107 ，所以 对于复化梯形求积公式有
3.复化 Gauss-Legendre I 型求积公式
将区间[a,b]划分为n等分，分点为xk a kh (k 0,1, 2, …,n;h 子区间[xk , xk 1 ]上采用2点Gauss-Legendre求积公式 (xk 1 xk ) x +x x -x 1 x +x x -x 1 [ f ( k 1 k - k 1 k )+f ( k 1 k + k 1 k )], 6 2 2 2 2 3 3 在[a,b]区间上的复化积分公式为 ba )，在每个 n
1 2
1
解:
一、复化求积公式基本介绍
1.复化梯形求积公式
将区间[a,b]划分为n等分，分点xk a kh( h 区间[ xk , xk 1 ](k 0,1, 2, …，n 1)上采用梯形式 I f ( x) dx
b a k 0 n 1 xk 1 xk b a
x 0
8
max( f 2(4) ) f 2(4)
1 2
x 0
96
又因为题目要求绝对误差限为 107 ，所以 对于复化梯形求积公式有
Rn ( f ) (b a) 2 '' (1 0)h 2 1 h f ( ) 8 = 10-7 12 12 2
所以
n 3651.5
ba , k 0,1, 2, …，n), 在每个子 n (b a ) f ( x) dx [ f ( a) f (b)], 则得 2
f ( x) dx
h [ f ( xk ) f ( xk 1)] Rn ( f ) 2 k 0
n 1
记
n 1 h n 1 h [ f ( xk ) f ( xk 1 )]= [ f ( a) 2 f ( xk ) f (b)] 2 k 0 2 k 1 称上式为复化梯形公式，其余项可由式
Байду номын сангаас
因此 取节点数 n=29 步长 h
ba 1 n 29
对于复化 Gauss-Legendre 型求积公式有
(b a)h 4 (4) (1 0)h 4 1 Rn ( f ) f ( ) 96 = 10-7 4320 4320 2
所以
n 25.8
因此 取节点数 n=26 步长 h
1 2 (3) 3x dx ln 3 0
3
2
1 1 1 d x (2) 4 dx 2 2 0 x 1 1 x
(4) e xe x dx
2 1
2
实验要求 : (1) 若用复化梯形公式、 复化 Simpson 公式和复化 Gauss-Legendre 型 公式做计算，要求绝对误差限为 107 ，分别利用它们的余项 对每种算法做出步长的事前估计。 (2) 分别用复化梯形公式、 复化 Simpson 公式和复化 Gauss-Legendre 型公式作计算。 (3) 将计算结果与精确解做比较，并比较各种算法的计算量。
Rn ( f ) (b a) 2 '' (1 0)h2 2 1 h f ( ) 4e = 10-7 12 12 2
所以
h ba 1 1.4248 10-4 n n n 7018.53
因此 取节点数 n=7019 步长 h 对于复化 Simpson 求积公式有
xk 1
xk
f ( x) dx
b
a
f ( x) dx
h n 1 h h [ f (x 1 )+f ( x 1 )] k k 2 k 0 2 3 2 3 2 2
上式称为复化Gauss-Legendre I 型求积公式。 ba 于是当f ( x) C 4 [a, b]，h 时，复化Gauss-Legendre I 型求积公式的余项表达式为 n (b a )h 4 (4) Rn ( f )= f ( ), [a, b] 4320