实验二 定积分的近似计算
定积分近似计算方法
定积分的近似计算方法摘要 本文主要讨论了一元函数常见的数值积分方法,例如插值型求积公式、龙贝格求积公式、高斯求积公式等近似计算方法,在用这些方法计算定积分时,会产生一些误差,为了减少误差, 可以利用复化求积公式、复化高斯公式等.本文围绕这些方法,系统介绍它们的计算公式以及截断误差,并用例题分析它们产生误差的大小、计算量等.关键词 插值型积分 龙贝格积分 高斯积分 误差分析 近似计算1引言在计算定积分的值()b aI f x dx =⎰时,常常根据微积分学基本定理求出)(x f 的一个原函数)(x F ,再用牛顿-莱布尼茨公式求的积分,()()()baI f x dx F b F a ==-⎰.但在实际应用中,这种方法只限于解决一小部分定积分的求值问题.当函数没有具体表达式,只是一些实验测得数据形成的表格或图形或者是()F x 无法用初等函数表示,例如,2bx ae dx ⎰,2sin ba x dx ⎰等等,这就需要我们用一些近似方法求的积分值.与数值积分一样,把积分区间细分,在每个小区间上,找到简单函数)(x ϕ来近似代替()f x ,且()bax dx ϕ⎰的值容易求的.这样就把计算复杂的()baf x dx ⎰转化为求简单的积分值()bax dx ϕ⎰.因此,定积分的近似计算实质上是就是被积函数的近似计算问题.2常见数值方法 2.1牛顿-科茨数值方法牛顿-科茨求积公式是求积节点等距离分布的插值型求积公式.利用插值多项式来构造数值积分公式是最常用、最基本的方法,具体做法是:给定区间[,]a b 上一组节点01...n a x x x b =<<<=,以及节点处函数()(0,1,2,i f x i n =,作()f x 的n 次拉格朗日多项式()()()nn i i i x f x l x ϕ==∑,其中 011011()()()()()()()()()i i n i i i i i i i n x x L x x x x L x x l x x x L x x x x L x x -+-+----=----,将插值公式(1)1()()()()(1)!n n n f f x x x n ξϕω++=++. 其中 1012()()()()()n n x xx x xx x L x x ω+=----,[,]a b ξ∈,依赖于变量x , 上式积分得(1)1()()()()(1)!n bb bn n aa af f x dx x dx x dx n ξϕω++=++⎰⎰⎰(1)(1)0()()()()(1)!n nb biiin aai f f x l x dx x dx n ξω++==++∑⎰⎰(1)(1)0()()()()(1)!n nb bi i n aai f f x b l x dx x dxn ξω++==++∑⎰⎰若记 (),(0,1,2,bi ia A l x dx i ==⎰….. )n (1)(1)1()[]()(1)!n bn af R f x dxn ξω++=+⎰, (2)则有()()[]nbi i ai f x dx A f x R f ==+∑⎰(3)称式(3)为插值求型公式,其中(0,1,2,i A i =…. )n 与()f x 无关,叫求积系数, i x 为求积节点,[]R f 为求积公式余项,其中求积系数由(1)决定.2.1.1梯形求积公式1梯形公式当插值节点01,x x 分别选取区间端点,a b 时,由式(3)分别求出求积系数10012bb aa x x xb b aA dx dx x x a b ---===--⎰⎰,01102bb aa x x x ab a A dx dx x x b a ---===--⎰⎰.从而的求积公式()[()()]2bab af x dx f a f b -≈+⎰. (4) 称求积公式(4)为梯形求积公式,简称梯形公式.2梯形公式截断误差: 3*()[](),12b a R f f ξ-''=- *[,]a b ξ∈. (5) 3梯形求积公式的代数精度:1 当()1f x =时,式(5)中 1(1)2bab adx b a x b a -=-=+=-⎰. 精确成立.2.1.2 辛普森求积公式1辛普森求积公式当选取节点为012,,2a bx a x x b +===时,由式(1)求下列求积系数 1200102()()()()2()()6()()2b b a a a b x x b x x x x b a A dx dx a b x x x x a a b +-----===+----⎰⎰,0211002()()()()2()()()3()()22bb aa x x x x x a xb b a A dx dx a b a b x x x x a b -----===++----⎰⎰.0122021()()()()2()()6()()22b b a a a bx a x x x x x b a A dx dx a b a b x x x x a b +-----===++----⎰⎰ .从而求积公式()[()4()()]62bab a a bf x dx f a f f b -+≈++⎰. (6)称式(6)为抛物线积分公式或辛普森积分公式.2抛物线求积公式误差估计定理1.若()f x 在[,]a b 上有四阶连续导数,则抛物线公式(6)的余项为:5(4)**()[](),[,]2880b a R f f a b ξξ--=∈. (7) 3抛物线公式的代数精度为3.易验证,当23()1,,,f x x x x =时,式(6)精确成立,而当4()f x x =时,式(6)不能精确成立.2.1.3 牛顿-科茨公式1牛顿-科茨公式在等距离节点i x a ih =+下,其中(0,1,2b ah i n-==…. )n .作为变量替换x a th =+,那么由求积公式(1),得系数:10(1)(1)(1)()!(1)(1)!ni n t t t i t i t n A h dt i n ---+---==--⎰10(1)(1)...(1)(1)...()(0,1,2,...)!(1)!n nb a t t t i t i t n dt i n n i n -----+---=-⎰ (8)则 ()()n i iA b a C =- (9)于是差值求积公式为:()0()()()[]nbn i i ai f x dx b a C f x R f ==-+∑⎰(10)称公式(10)为牛顿-科茨求积公式,其中()n iC 称为科茨系数.显然,科茨系数与被积函数()f x 及积分区间[,]a b 无关,它指依赖于n ,且为多项式积分.因此,只要给出n ,就能看出i A ,并写出相应地牛顿-科茨公式.2牛顿-科茨公式的截断误差与代数精度.当1n =与2n =情况分析牛顿-科茨公式的截断误差为(1)()[]()()()(1)!n b b bn aaaf R f f x dx x dx x dxn ξϕω+=-=+⎰⎰⎰牛顿-科茨公式的截断误差还可以写成(2)*1()[]()((2)!n bn a f R f x dx n n ξω++=+⎰为偶数)(1)*1()[]()(1)!n bn af R f x dx n ξω++=+⎰ (n 为奇数) (11) 其中*[,]a b ξ∈,且不依赖于x ,101()()()...()n n x x x x x x x ω+=---,对()f x 为任何并不超过n 次多项式,均有(1)()0n fx +≡,因而[]0R f ≡,即0()()nbi i ai f x dx A f x ==∑⎰精确成立,也就是说,牛顿-科茨公式的代数精度至少为n ,牛顿-科茨公式在n 为偶数时,至少具有1n +次代数精度,在n 为奇数情况时,至少具有n 次代数精度.2.1.4复化梯形求积公式将区间[,]a b 等分,节点为i x a ih =+ (步长b ah n-=),0,1,2...,i n =)在每个小区间1[,]i i x x -上采用梯形公式(4)得11111()()[(()()]2ii nnbx i i i i ax i i x x f x dx f x dx f x f x ---==-=≈+=∑∑⎰⎰11[()()]2ni i i hf x f x +=+=∑11[()2()()]2n i n i hf a f x f b T -=++=∑ (12)称式(12)为复化梯形公式. 复化梯形公式余项为()2()()()12i n b a R f h f η-''=-(13) 2.1.5复化辛普森求积公式在每个小区间],[1+i i x x 上,辛普森公式(6)得11102()[()4()()]6n bi i ai i hf x dx f x f x f x -++==++∑⎰(14)111012[()4()2((6)]6n n i i i i hf a f x f x f --+===+++∑∑记 )]()(2)(4)([6111021b f x f x f a f hS n i i n i i n +++=∑∑-=-=+ (15)式中,21+i x为],[1+i i x x 的中点,即h x x i i 2121+=+.式(15)称为复化辛普森公式,其余项为∑-=-=-=10)4(4)()2(180)()(n i i n n f h h S f I f R η, 1(,).i i i x x η+∈ 故 ),(),()2(180)(R )4(4b a f h a b f n ∈--=ηη (16) 为复化辛普森的截断误差. 2.1.6复化科茨求积公式将区间[,]a b n 等分, 4n m =,m 为正整数,在每个子区间444[,]k k x x -上用科茨求积公式得到复化求积公式:412()[7()7()32()45mbk ak hf x dx f a f b f x -≈++∑⎰14241411112()32()14()mmm k k k N k k k f xf x f x C ---===+++=∑∑∑ (17)其中 4b a b ah n m--==, k x a kh =+ 其截断误差为6(6)2()[,](),()945n b a R f C h f a b ηη-=-<. 2.1.7 变步长复化求积方法复化求积公式虽然计算简单,也达到了提高精度的目的,但为了满足精度要求必须顾及误差,利用误差公式往往很困难,因为误差表达式中含有未知函数的导数,而估计各阶导数的最大值不太容易.我们可以采取把积分的区间[,]a b 细分的办法,在计算积分时将步长逐步折半,利用前后两次结果进行误差估计,如此继续,直到相邻两次结果相差不大,取最小的步长算出的结果为积分值,这种方法称为变步长积分法.以复化梯形公式为例,把区间[,]a b 分成n 等分,设复化梯形公式的近似值为n T ,原积分值为I ,由复化梯形公式误差公式(14)知:2"11()()()n b a b a I T f a b N N ηη--=-<<再把区间[,]a b 分成2n 等分,得近似值2n T ,则2222()()()122k b a b a I T f a b nηη--''=-<< 假定()f x ''在[,]a b 上变化不大,既有12()()f f ηη''''≈. 由上式得 .24kkI T I T -≈-于是 222211()()341n n n n n n I T T T T T T ≈+-=+-- (18) 式(18)表明若用2n T 作为I 的近似值,其截断误差约为2()3n n T T - (19)2.2 龙贝格求积公式龙贝格积分法的基本思想是采用复化梯形求积方法不断折半步长过程中,在积分结果中加入时候误差估计值进行补偿,使积分计算的收敛性加速,就可以加工出,,,...n n n S C R 精度较高的积分结果.由式(19), 2n T 的误差大致为23n nT T -,因此,可用这个误差值作为2n T 的一种补偿,加到2n T 上,则可得到积分准确值I ,比2n T 的更好近似值~T .222141()333n n n n nT T T T T T =+-=-2221(2)21n n T T =-- (20)式(20)左端1n =时 记122121141()333S T T T T T =+-=- 112()()332a b T b a f +=+- [()4()()]62b a a b f a f f b -+=++恰好为[,]a b 上应用辛普生公式(16)的结果.在每个小区间应用辛普生公式:11[()2()()]2n n k k hT f a f x f b -==++∑121()112[()2()()2()]4n n n k k k k hT f a f x f b f x --===+++∑∑代入式(20)的左端得11111[()2()()2()32n nk k k k h f a f x f b f x -==+++--∑∑ 11[()2()()]2n k k h f a f x f b -++∑11111[()4()2()()]62n n k k k k f a f x f x f b -===+-++∑∑nS =从而复化辛普森公式与复化梯形公式公式有以下关系式2441n nn T T S -=- (21)类似也可以推证,在辛普森序列基础上,利用以下关系式22242161151541n n n n n S S C S S -=-=- (22)可以造出收敛速度更快的科茨序列12,...,...n C C C 将此推行下去,在科茨序列基础上,通过243431n nn C C R -=- (23)构造出收敛速度比科茨序列更快的龙贝格序列12,,......n R R R .以上这种通过逐步构造龙贝格序列的积分近似值法就称为龙贝格积分法.2.3高斯求积公式由定理()()()baf x F b F a =-⎰知,插值型求积公式的代数精度与求积节点的个数有关,具有1n +个节点的插值型求积公式至少具有n 次代数精度.不仅如此,代数精度与节点的选取有关,在构造牛顿-科茨求积公式时,为了简化处理过程,限定用等分节点作为求积节点,这样做,虽然公式确实得到简化,但同时也限制了公式的代数精度. 设积分,1,1=-=b a 本段讨论如下求积公式11()()ni i i f x A f x -==∑⎰(24)对任意积分区间[,]a b ,通过变 22ba t ab x ++-= 可以转换到区间]1,1[-上,这时11()()222bab a b a a bf x dx f t dt ---+=+⎰⎰ 此时,求积公式写为0()()222n bii ai b a a b b af x dx A f t =-+-=+∑⎰若一组节点]1,1[.....,10-∈n x x x 使插值型求积公式(24)具有21n +次代数精度,则称此组节点为高斯点,并称相应求积公式(24)为高斯求积公式.2.3.1 高斯求积公式的余项(2)2()[]()()()(22)!n nbb k k aa k f R f f x dx A f x x dx n ηω+==-=+∑⎰⎰ 其中 01()()()...(),[,]n x x x x x x x ab ωη=---∈,且不依赖于x .2.3.2 复化高斯求积公式复化高斯求积公式的基本思想是:将积分区间[,]a b 分成n个等长小区间1[,](1,...)i i t t i m -=,然后在低阶(2n =)高斯求积公式算出近似值,最后将他们相加的积分()baf t dt ⎰的近似值m G ,即11111111()()[]222ii mmbt i i i i i i at i i t t t t t t f t dt f t dt dt -----==-+-==+∑∑⎰⎰⎰1111[()]222m i h ha i h x dx-==+-+∑⎰101[()]222m n j j mi j h hA f a i h x G ==≈+-+≈∑∑ (25)其中mab h -=,j A 与(0,1,2,...,)j t j n =可由书中表中查出. 3 应用3.1插值型积分的应用例1 用牛顿-科茨公式(1,2,4n =)计算积分12211I x =+⎰. 解 1n =时2210112[]0.4512101()2I -≈+=++2n =时22211112[4]0.463725116101()1()42I -≈++=+++4n =时2222111112[7321232]0.46363311390101()1()1()848I =++++≈++++例2 利用复化梯形求积公式计算积分 12211I dx x =+⎰解 设211)(xx f +=,分点个数为n =1,2,4,5时,求出相应积分n T , 111[(()())],21,2(),.n n i i i i i T f a f b f h b a h n n f x f x a ih ih -=⎧=++⎪⎪-⎪==⎨⎪=⎪⎪=+=⎩∑列表如下:n =1的计算结果见表1-1所列 n h0x 1x 0f1f1T10.50.00.51.0 0.8 0.45n =2的表格如下 n hx1x2xf1f2f2T20.250.00 0.25 0.50 1.00 0.941765 0.80 0.460294n =4时计算结果如下表 n h 0x1x2x3x4x40.1250.00 0.125 0.25 0.375 0.50f1f2f3f4f4T1.00 0.9846154 0.9411765 0.876712 0.80 0.462813n = 5时计算结果如下 n hx1x2x3x4x5x50.10.0 0.1 0.2 0.3 0.40.5f1f2f3f4f5f5T1.0 0.990099 0.9615385 0.91743 0.862069 0.80.463114例3 利用复化求积公式120x e dx ⎰,问积分区间为多少等分才能得证有5位有效数字?解 由式(14)知322()[],()()1212n b a b a R f h f n f n n--''''=-=- 有1(),(),2x x f x e f x e b a ''==-=,当]21,0[∈x 时,在12|()|f x e ''≤,所以122|[]|96n eR f n≤ 由于120x e dx ⎰的准确值具有一位整数,所以要使近似值具有5位有效数字,n 必须满足4242211048,102196⨯≥⨯≤-e n n e 或 取对数有 19=n .即将区间]21,0[19等分可满足给定的精度要求.例4 利用复化抛物线求积公式计算 120211I dx x =+⎰. 解 设11)(2+=x x f ,取m =1,2, 3时,公式()⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧++=+=====-=+++=+---=-=+∑∑.)12(,2),(),(),(,,242[31221212221111,1222h i a x ih a x x f f x f f b f f a f f m a b n f f f f h S i i i i i i b a m i m i i b a m当m =1,2,3时结果如下表所示 当m =1时m h(0.0)f )25.0(f )5.0(f2S1 0.25 1.0 0.9411765 0.80 0.463725当m =2时mh(0.0)f(0.125)f (0.025)f (0.35)f )5.0(f4S20.125 1.0 0.9846154 0.9411765 0.8767123 0.80 0.463653当m =3时mh(0.0)f(0.08333)f (0.16667)f (0.35)f(0.33333)f (0.14166667)f )5.0(f4S30.83331.00.99310340.9729730.9411760.90.852070.80.4636例5 用复化梯形公式,辛普森公式和科茨公式计算积分10sin xdx x ⎰的近似值.解按精度要求确定]1,0[分多少等分,即确定步长,要使6441021)1(28801|],[|-⨯≤≤M m S f R n ,只需.4642880102M m ⨯≥令10sin ()cos xf x txdt x==⎰,则1()0sin ()()(cos )k kk k k d xd fx tx dt dx x dx==⎰ 1cos().2k t tx kdt π=+⎰dt ktx t x f k k |)2cos(|max )(|max 10)(π+≤⎰11.1k t d t t≤=+⎰)10(≤≤x (4)1max |()| 5.f x ≤所以只要,9.13831288010264=⨯⨯≥-m 取m =4即可, 当4n =时,在每个子区间上用式(25),或(14),或(17),结果.9460829.0,9460833.0,9456911.0888===C S T3.2 龙贝格积分公式应用例6 用龙贝格算法计算积分1241I dx x=+⎰的近似值,要求误差小于510-. 解 .3,0,14)(2==+=b a x x f 步骤如下:2)1(,4)0()1(==f f 得.3)]1()0([211=+=f f T )2(计算,1.3)]21([21,516)21(12=+==f T T f 由此得301333334121=-=T T S . (3)算出),(43),41(f f 从而,3013118)]43()41([412124=++=f f T T,14157.334242=-=T T S .30142121516121=-=S S C(4)计算),87(),85(),83(),81(f f f f 从而得到:13899.3)]87()85()83()81([812148=++++=f f f f T T ,,14159.334482=-=T T S ,14059.31516242=-=S S C.1458.36364121=-=C C R (5)再计算),1615(),1613(),1611(),169(),167(),165(),163(),161(f f f f f f f f 从而得到: 14094.316=T30141598=S ,,14159.3,14159.324==R C 51210||-≤-R R , 所以12043.14159.1dx x ≈+⎰3.3高斯求积公式的应用例7 用两点复化高斯求积公式计算10,x I e dx =⎰要求允许误差.106-=ε解 在本算法中取21=+n 时,,110==A A 其中;,)(mab h e x f x-== =++--=∑=)22(2201j jj b a x a b f A a b G.87189637800.1][21)32121()32121(=++-eem =2时, h =21, ]4121)21([4120202j i j j x i f A G +⨯-=∑∑==.57182571650.1)(41341333413341333413=+++=++--eeee m =3时, h =31. .37182769352.1]631)21([6130203=+⨯-=∑∑==j i j j x i f A G.101027.71||||56323--<⨯≈+-G G G3.4 几种方法的比较分析例8 计算积分211ln 2dx x =⎰,精确到0.001.(1)利用矩形公式计算, 因为对于x x f 1)(=,有320()2f x x''<=<(如果1<x <2),所以按照公式0)2(S =+-dx ba xb a . 0<n R <2112n . 如果取n =10,则我们公式的余项的余数得31010.84101200R -<<⨯,我们还必须加进由于在计算函数值实行四舍五入所产生的误差的界限相差于0.16⨯310-,为了这个目的只要计算1x的值到四位小数精确到0.00005就够了.我们有1232527292132152172192 1.051.151.251.351.551.651.751.851.95x x x x x x x x x =========5128.05405.05714.06061.06897.07407.08.08696.09524.02192172152132927252321=========y y y y y y y y y和6.928469284.0109284.6= (2) 按照梯形公式作同样的计算,在这种情况下,作公式 210,||6n n R R n<<在这儿也试一试取n =10,虽然此时仅可以证3107.16001||-⨯<<n R ,纵坐标是9.18.17.16.15.14.13.12.11.1987654321=========x x x x x x x x x 5263.05556.05882.06250.06667.07143.07692.08333.09091.0987654321=========y y y y y y y y y和1877.669377.01877.621500101=+)( (3) 用辛普森公式做同样的计算作公式 .0))(()2(180)()4(45<≤≤⨯--=n n R b a f n a b R ξξ 并且n =5时有55104.1||-⨯<R .实行计算到五位数字,精确到0.0000058.16.14.12.14321====x x x x 45636.555556.062500.071429.083333.04321和====y y y y 9.17.15.13.11.12927252321=====x x x x x83820.1352632.058824.066667.076923.090909.029********和=====y y y y y.20.150==x x 50000.150000.060000.150和==y y6931525.083820.345636.550000.1301=++)(. 由此可见,用辛普森公式计算得到的值误差最小,计算量相对一般;而用矩形公式计算得到的值误差较大,计算量也比较大;用梯形公式计算的值误差比用矩形公式得到的值要误差小,计算量也是如此.所以我们计算定积分时用辛普森公式往往得到的值误差小,而对没有要求误差大小的,则可以选择辛普森或者是梯形公式,因为这两种方法计算量相对较小.结 束 语本文只讨论了一些一维数值积分方法及其它们的应用,误差分析等有关内容.其中最常用的方法是插值型积分以及复化方法、龙贝格积分方法和高斯积分方法,并讨论了相关求积方法的代数精度和误差分析,并给出了一些例题,分析各种方法的近似值,得出误差分析最小的近似方法.由于篇幅有限,对于高维数值积分方法本文便不再讨论.参考文献[1] 华东师范大学数学系,数学分析(第一版)[M],北京:高等教育出版社,2001. [2] 李庆阳,关治,白峰杉,数值计算原理(第二版)[M],北京: 清华大学出版社, 2008. [3] 肖筱南,现代数值计算方法(第一版)[M],北京: 北京大学出版社, 1999.[4] 菲赫金格尔茨,微积分学教程(第三版)[M],北京: 高等教育出版社, 2005. [5] 裴礼文,数学分析中的典型问题与方法(第一版)[M] ,北京: 北京大学出版社,2004. [6] 李桂成,计算方法(第三版)[M],北京: 高等教育出版社,2010.[7] Yin Y uezhu ,Yang Zhonglian.Calculating Skillfully the Curve Integral and Surface Integral Type 2 bySymmetry, SCIENCE & TECHNOLOGY INFORMATION ,2008(30)The Approximate Numerical Method of the Definite IntegralAbstract This paper mainly discusses common numerical methods of unary function, such as approximate calculation method of interpolation integral, Lebesgue integral and Gauss integration. With these methods in calculating the integral, it will produce some error. In order to reduce the error, we can use after the formula for product and after the Gauss formula. This paper focus on these methods introducing formula of introduction and truncation errors .In addition they can provide examples to analysis size of the error and computation.Keywords interpolation integral Lebesgue integral Gauss integral error analysis approximate computation。
《数学实验》实验报告——定积分的近似求解
3
2 梯形法程序如下: f=input('请输入被积函数 f(x)='); qujian=input('请输入积分区间[a,b]='); n=input('请输入子区间个数 n='); s=0; for i=1:n-1 x=qujian(1)+(qujian(2)-qujian(1))/n*i; y=eval(f); s=s+y; end x=qujian(1); y=eval(f); s=s+y/2; x=qujian(2); y=eval(f); s=s+y/2; disp('定积分的近似值是:'); s=s*(qujian(2)-qujian(1))/n
《数学实验》实验报告
班级 试验 内容 **** 学号 **** 姓名 试验 类别
自选试验
****
成绩 试验 时间 2011 年 5 月 20 日—22 日
定积分的近似求解
试验问题:
用梯形法与抛物法,通过 MATLAB,计算 x 2 dx 的近似值,取 n=10,比较结果的差异,研究
0 1
定积分的两种近似计算方法。
1 1 1 2 ph 3 6rh h(2 ph 2 6r ) h( y 0 4 y1 y 2 ) 3 3 3 。 ba n ,则上面所求的 S 等于区间 [ x0 , x2 ] 上以抛物线为曲边的曲边梯形的面积。同理可
取
以得到区间 [ xi 1 , xi 1 ] 上以抛物线为曲边的曲边梯形的面积:
试验目的:
通过分别用梯形法与抛物线法计算定积分的近似值, 进而熟练掌握运用 MATLAB 来解决 定积分的近似求解,体会 MATLAB 的强大功能。
matlab实验报告--定积分的近似计算
○2 使用函数 quad()
quad('sin(x)./x',0,inf) 【调试结果】 ans =
NaN
○3 程序法
%矩阵法
format long
n=inf;a=0;b=inf;
syms x fx
fx=sin(x)./x;
i=1:n;
xj=a+(i-1)*(b-a)/n; xi=a+i*(b-a)/n;
%左点
xi=a+i*(b-a)/n;
%右点
xij=(xi+xj)/2;
fxj=subs(fx,'x',xj);
%左点值
fxi=subs(fx,'x',xi);
%右点值
fxij=subs(fx,'x',xij);
%中点值
f=(fxj+4*fxij+fxi)*(b-a)/(6*n);
inum=sum(f)
xj=a+(i-1)*(b-a)/n; xi=a+i*(b-a)/n; xk=(xi+xj)/2;
%左点 %右点 %中点
fxj=subs(fx,'x',xj); fxi=subs(fx,'x',xi); fxk=subs(fx,'x',xk); inum=inum+(fxj+4*fxk+fxi)*(b-a)/(6*n); end inum integrate=int(fx,1,2); integrate=double(integrate); fprintf('The relative error between inum and real-value is about:%g/n/n',... abs((inum-integrate)/integrate)) 【调试结果】 inum =
定积分近似计算方法
定积分的近似计算方法摘要 本文主要讨论了一元函数常见的数值积分方法,例如插值型求积公式、龙贝格求积公式、高斯求积公式等近似计算方法,在用这些方法计算定积分时,会产生一些误差,为了减少误差, 可以利用复化求积公式、复化高斯公式等.本文围绕这些方法,系统介绍它们的计算公式以及截断误差,并用例题分析它们产生误差的大小、计算量等.关键词 插值型积分 龙贝格积分 高斯积分 误差分析 近似计算1引言在计算定积分的值()b aI f x dx =⎰时,常常根据微积分学基本定理求出)(x f 的一个原函数)(x F ,再用牛顿-莱布尼茨公式求的积分,()()()baI f x dx F b F a ==-⎰.但在实际应用中,这种方法只限于解决一小部分定积分的求值问题.当函数没有具体表达式,只是一些实验测得数据形成的表格或图形或者是()F x 无法用初等函数表示,例如,2bx ae dx ⎰,2sin ba x dx ⎰等等,这就需要我们用一些近似方法求的积分值.与数值积分一样,把积分区间细分,在每个小区间上,找到简单函数)(x ϕ来近似代替()f x ,且()b a x dx ϕ⎰的值容易求的.这样就把计算复杂的()ba f x dx ⎰转化为求简单的积分值()bax dx ϕ⎰.因此,定积分的近似计算实质上是就是被积函数的近似计算问题.2常见数值方法 2.1牛顿-科茨数值方法牛顿-科茨求积公式是求积节点等距离分布的插值型求积公式.利用插值多项式来构造数值积分公式是最常用、最基本的方法,具体做法是:给定区间[,]a b 上一组节点01...n a x x x b =<<<=,以及节点处函数()(0,1,2,i f x i n =,作()f x 的n 次拉格朗日多项式()()()nn i i i x f x l x ϕ==∑,其中 011011()()()()()()()()()i i n i i i i i i i n x x L x x x x L x x l x x x L x x x x L x x -+-+----=----,将插值公式(1)1()()()()(1)!n n n f f x x x n ξϕω++=++.其中1012()()()()()n n x x x x x x x L x x ω+=----,[,]a b ξ∈,依赖于变量x , 上式积分得(1)1()()()()(1)!n bb bn n aa af f x dx x dx x dx n ξϕω++=++⎰⎰⎰(1)(1)0()()()()(1)!n nb biiin aai f f x l x dx x dx n ξω++==++∑⎰⎰(1)(1)0()()()()(1)!n nbbi i n aai f f x b l x dx x dxn ξω++==++∑⎰⎰若记 (),(0,1,2,bi ia A l x dx i ==⎰….. )n (1)(1)1()[]()(1)!n bn af R f x dxn ξω++=+⎰, (2)则有()()[]nbi i ai f x dx A f x R f ==+∑⎰(3)称式(3)为插值求型公式,其中(0,1,2,i A i =…. )n 与()f x 无关,叫求积系数, i x 为求积节点,[]R f 为求积公式余项,其中求积系数由(1)决定.2.1.1梯形求积公式1梯形公式当插值节点01,x x 分别选取区间端点,a b 时,由式(3)分别求出求积系数10012bb aa x x xb b aA dx dx x x a b ---===--⎰⎰,01102bb aa x x x ab a A dx dx x x b a ---===--⎰⎰.从而的求积公式()[()()]2bab a f x dx f a f b -≈+⎰. (4) 称求积公式(4)为梯形求积公式,简称梯形公式.2梯形公式截断误差: 3*()[](),12b a R f f ξ-''=- *[,]a b ξ∈. (5) 3梯形求积公式的代数精度:1 当()1f x =时,式(5)中 1(1)2bab adx b a x b a -=-=+=-⎰. 精确成立.2.1.2 辛普森求积公式1辛普森求积公式当选取节点为012,,2a bx a x x b +===时,由式(1)求下列求积系数 1200102()()()()2()()6()()2b b a a a b x x b x x x x b a A dx dx a b x x x x a a b +-----===+----⎰⎰,0211002()()()()2()()()3()()22bb aa x x x x x a xb b a A dx dx a b a b x x x x a b -----===++----⎰⎰.0122021()()()()2()()6()()22b b a a a bx a x x x x x b a A dx dx a b a b x x x x a b +-----===++----⎰⎰ .从而求积公式()[()4()()]62bab a a bf x dx f a f f b -+≈++⎰. (6)称式(6)为抛物线积分公式或辛普森积分公式.2抛物线求积公式误差估计定理1.若()f x 在[,]a b 上有四阶连续导数,则抛物线公式(6)的余项为:5(4)**()[](),[,]2880b a R f f a b ξξ--=∈. (7) 3抛物线公式的代数精度为3.易验证,当23()1,,,f x x x x =时,式(6)精确成立,而当4()f x x =时,式(6)不能精确成立.2.1.3 牛顿-科茨公式1牛顿-科茨公式在等距离节点i x a ih =+下,其中(0,1,2b ah i n-==…. )n .作为变量替换x a th =+,那么由求积公式(1),得系数:10(1)(1)(1)()!(1)(1)!ni n t t t i t i t n A h dt i n ---+---==--⎰10(1)(1)...(1)(1)...()(0,1,2,...)!(1)!n nb a t t t i t i t n dt i n n i n -----+---=-⎰ (8)则 ()()n i i A b a C =- (9) 于是差值求积公式为:()0()()()[]nbn i i ai f x dx b a C f x R f ==-+∑⎰(10)称公式(10)为牛顿-科茨求积公式,其中()n iC 称为科茨系数.显然,科茨系数与被积函数()f x 及积分区间[,]a b 无关,它指依赖于n ,且为多项式积分.因此,只要给出n ,就能看出i A ,并写出相应地牛顿-科茨公式.2牛顿-科茨公式的截断误差与代数精度.当1n =与2n =情况分析牛顿-科茨公式的截断误差为(1)()[]()()()(1)!n b b bn aaaf R f f x dx x dx x dxn ξϕω+=-=+⎰⎰⎰牛顿-科茨公式的截断误差还可以写成(2)*1()[]()((2)!n bn a f R f x dx n n ξω++=+⎰为偶数)(1)*1()[]()(1)!n bn a f R f x dx n ξω++=+⎰ (n 为奇数) (11) 其中*[,]a b ξ∈,且不依赖于x ,101()()()...()n n x x x x x x x ω+=---,对()f x 为任何并不超过n 次多项式,均有(1)()0n fx +≡,因而[]0R f ≡,即0()()nbi i ai f x dx A f x ==∑⎰精确成立,也就是说,牛顿-科茨公式的代数精度至少为n ,牛顿-科茨公式在n 为偶数时,至少具有1n +次代数精度,在n 为奇数情况时,至少具有n 次代数精度.2.1.4复化梯形求积公式将区间[,]a b 等分,节点为i x a ih =+ (步长b ah n-=),0,1,2...,i n =)在每个小区间1[,]i i x x -上采用梯形公式(4)得11111()()[(()()]2ii nnbx i i i i ax i i x x f x dx f x dx f x f x ---==-=≈+=∑∑⎰⎰11[()()]2ni i i hf x f x +=+=∑11[()2()()]2n i n i hf a f x f b T -=++=∑ (12)称式(12)为复化梯形公式. 复化梯形公式余项为()2()()()12i n b a R f h f η-''=-(13) 2.1.5复化辛普森求积公式在每个小区间],[1+i i x x 上,辛普森公式(6)得11102()[()4()()]6n bi i ai i hf x dx f x f x f x -++==++∑⎰(14)111012[()4()2((6)]6n n i i i i hf a f x f x f --+===+++∑∑记 )]()(2)(4)([6111021b f x f x f a f hS n i i n i i n +++=∑∑-=-=+ (15)式中,21+i x为],[1+i i x x 的中点,即h x x i i 2121+=+.式(15)称为复化辛普森公式,其余项为∑-=-=-=10)4(4)()2(180)()(n i i n n f h h S f I f R η, 1(,).i i i x x η+∈故 ),(),()2(180)(R )4(4b a f h a b f n ∈--=ηη (16) 为复化辛普森的截断误差. 2.1.6复化科茨求积公式将区间[,]a b n 等分, 4n m =,m 为正整数,在每个子区间444[,]k k x x -上用科茨求积公式得到复化求积公式:412()[7()7()32()45mbk ak hf x dx f a f b f x -≈++∑⎰14241411112()32()14()mmm k k k N k k k f xf x f x C ---===+++=∑∑∑ (17)其中 4b a b a h n m--==, k x a kh =+ 其截断误差为6(6)2()[,](),()945n b a R f C h f a b ηη-=-<. 2.1.7 变步长复化求积方法复化求积公式虽然计算简单,也达到了提高精度的目的,但为了满足精度要求必须顾及误差,利用误差公式往往很困难,因为误差表达式中含有未知函数的导数,而估计各阶导数的最大值不太容易.我们可以采取把积分的区间[,]a b 细分的办法,在计算积分时将步长逐步折半,利用前后两次结果进行误差估计,如此继续,直到相邻两次结果相差不大,取最小的步长算出的结果为积分值,这种方法称为变步长积分法.以复化梯形公式为例,把区间[,]a b 分成n 等分,设复化梯形公式的近似值为n T ,原积分值为I ,由复化梯形公式误差公式(14)知:2"11()()()n b a b a I T f a b N N ηη--=-<<再把区间[,]a b 分成2n 等分,得近似值2n T ,则2222()()()122k b a b a I T f a b nηη--''=-<< 假定()f x ''在[,]a b 上变化不大,既有12()()f f ηη''''≈. 由上式得 .24kkI T I T -≈-于是 222211()()341n n n n n n I T T T T T T ≈+-=+-- (18) 式(18)表明若用2n T 作为I 的近似值,其截断误差约为2()3n n T T - (19)2.2 龙贝格求积公式龙贝格积分法的基本思想是采用复化梯形求积方法不断折半步长过程中,在积分结果中加入时候误差估计值进行补偿,使积分计算的收敛性加速,就可以加工出,,,...n n n S C R 精度较高的积分结果.由式(19), 2n T 的误差大致为23n nT T -,因此,可用这个误差值作为2n T 的一种补偿,加到2n T 上,则可得到积分准确值I ,比2n T 的更好近似值~T .222141()333n n n n nT T T T T T =+-=- 2221(2)21n n T T =-- (20)式(20)左端1n =时 记122121141()333S T T T T T =+-=- 112()()332a b T b a f +=+- [()4()()]62b a a b f a f f b -+=++恰好为[,]a b 上应用辛普生公式(16)的结果.在每个小区间应用辛普生公式:11[()2()()]2n n k k hT f a f x f b -==++∑121()112[()2()()2()]4n n n k k k k hT f a f x f b f x --===+++∑∑代入式(20)的左端得11111[()2()()2()32n nk k k k h f a f x f b f x -==+++--∑∑ 11[()2()()]2n k k h f a f x f b -++∑11111[()4()2()()]62n n k k k k f a f x f x f b -===+-++∑∑n S =从而复化辛普森公式与复化梯形公式公式有以下关系式2441n nn T T S -=- (21)类似也可以推证,在辛普森序列基础上,利用以下关系式22242161151541n n n n n S S C S S -=-=- (22)可以造出收敛速度更快的科茨序列12,...,...n C C C 将此推行下去,在科茨序列基础上,通过243431n nn C C R -=- (23)构造出收敛速度比科茨序列更快的龙贝格序列12,,......n R R R .以上这种通过逐步构造龙贝格序列的积分近似值法就称为龙贝格积分法.2.3高斯求积公式由定理()()()baf x F b F a =-⎰知,插值型求积公式的代数精度与求积节点的个数有关,具有1n +个节点的插值型求积公式至少具有n 次代数精度.不仅如此,代数精度与节点的选取有关,在构造牛顿-科茨求积公式时,为了简化处理过程,限定用等分节点作为求积节点,这样做,虽然公式确实得到简化,但同时也限制了公式的代数精度. 设积分,1,1=-=b a 本段讨论如下求积公式11()()ni i i f x A f x -==∑⎰(24)对任意积分区间[,]a b ,通过变 22ba t ab x ++-= 可以转换到区间]1,1[-上,这时11()()222bab a b a a bf x dx f t dt ---+=+⎰⎰ 此时,求积公式写为0()()222n bii ai b a a b b af x dx A f t =-+-=+∑⎰若一组节点]1,1[.....,10-∈n x x x 使插值型求积公式(24)具有21n +次代数精度,则称此组节点为高斯点,并称相应求积公式(24)为高斯求积公式.2.3.1 高斯求积公式的余项(2)2()[]()()()(22)!n nbb k k aa k f R f f x dx A f x x dx n ηω+==-=+∑⎰⎰ 其中01()()()...(),[,]n x x x x x x x a b ωη=---∈,且不依赖于x .2.3.2 复化高斯求积公式复化高斯求积公式的基本思想是:将积分区间[,]a b 分成n 个等长小区间1[,](1,...)i i t t i m -=,然后在低阶(2n =)高斯求积公式算出近似值,最后将他们相加的积分()baf t dt ⎰的近似值m G ,即11111111()()[]222ii mmbt i i i i i i at i i t t t t t t f t dt f t dt dt -----==-+-==+∑∑⎰⎰⎰1111[()]222m i h ha i h x dx-==+-+∑⎰101[()]222m n j j mi j h hA f a i h x G ==≈+-+≈∑∑ (25)其中mab h -=,j A 与(0,1,2,...,)j t j n =可由书中表中查出. 3 应用3.1插值型积分的应用例1 用牛顿-科茨公式(1,2,4n =)计算积分12211I x =+⎰. 解 1n =时2210112[]0.4512101()2I -≈+=++2n =时22211112[4]0.463725116101()1()42I -≈++=+++4n =时2222111112[7321232]0.46363311390101()1()1()848I =++++≈++++例2 利用复化梯形求积公式计算积分 12211I dx x =+⎰解 设211)(xx f +=,分点个数为n =1,2,4,5时,求出相应积分n T , 111[(()())],21,2(),.n n i i i i i T f a f b f h b a h n n f x f x a ih ih -=⎧=++⎪⎪-⎪==⎨⎪=⎪⎪=+=⎩∑列表如下:n =1的计算结果见表1-1所列 n h0x 1x 0f1f1T10.50.00.51.00.80.45n =2的表格如下 n h0x1x2x0f1f 2f 2T20.250.000.250.501.000.941765 0.800.460294n =4时计算结果如下表 n h 0x1x2x3x4x40.1250.000.1250.250.3750.500f1f2f3f4f4T1.000.98461540.94117650.8767120.800.462813n = 5时计算结果如下 n h0x1x2x3x4x5x50.10.00.10.20.30.40.50f1f2f3f4f5f5T1.00.9900990.96153850.917430.8620690.80.463114例3 利用复化求积公式120x e dx ⎰,问积分区间为多少等分才能得证有5位有效数字?解 由式(14)知322()[],()()1212n b a b a R f h f n f n n--''''=-=- 有1(),(),2x xf x e f x e b a ''==-=,当]21,0[∈x 时,在12|()|f x e ''≤,所以122|[]|96n eR f n≤ 由于120x e dx ⎰的准确值具有一位整数,所以要使近似值具有5位有效数字,n 必须满足4242211048,102196⨯≥⨯≤-e n n e 或 取对数有 19=n .即将区间]21,0[19等分可满足给定的精度要求.例4 利用复化抛物线求积公式计算 120211I dx x =+⎰. 解 设11)(2+=x x f ,取m =1,2, 3时,公式()⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧++=+=====-=+++=+---=-=+∑∑.)12(,2),(),(),(,,242[31221212221111,1222h i a x ih a x x f f x f f b f f a f f m a b n f f f f h S i i i i i i b a m i m i i b a m当m =1,2,3时结果如下表所示 当m =1时m h(0.0)f )25.0(f )5.0(f2S1 0.25 1.0 0.9411765 0.80 0.463725当m =2时mh(0.0)f(0.125)f (0.025)f (0.35)f )5.0(f4S20.125 1.0 0.9846154 0.9411765 0.8767123 0.80 0.463653当m =3时mh(0.0)f(0.08333)f (0.16667)f (0.35)f(0.33333)f (0.14166667)f )5.0(f4S30.83331.00.99310340.9729730.9411760.90.852070.80.4636例5 用复化梯形公式,辛普森公式和科茨公式计算积分10sin xdx x ⎰的近似值.解按精度要求确定]1,0[分多少等分,即确定步长,要使6441021)1(28801|],[|-⨯≤≤M m S f R n ,只需.4642880102M m ⨯≥令10sin ()cos xf x txdt x==⎰, 则1()0sin ()()(cos )k kk k k d xd fx tx dt dx xdx ==⎰1cos().2k t tx kdt π=+⎰dt ktx t x f k k |)2cos(|max )(|max 10)(π+≤⎰11.1k t d t t≤=+⎰)10(≤≤x (4)1max |()| 5.f x ≤所以只要,9.13831288010264=⨯⨯≥-m 取m =4即可, 当4n =时,在每个子区间上用式(25),或(14),或(17),结果.9460829.0,9460833.0,9456911.0888===C S T3.2 龙贝格积分公式应用例6 用龙贝格算法计算积分1241I dx x=+⎰的近似值,要求误差小于510-. 解 .3,0,14)(2==+=b a x x f 步骤如下:2)1(,4)0()1(==f f 得.3)]1()0([211=+=f f T )2(计算,1.3)]21([21,516)21(12=+==f T T f 由此得 301333334121=-=T T S . (3)算出),(43),41(f f 从而,3013118)]43()41([412124=++=f f T T,14157.334242=-=T T S.30142121516121=-=S S C (4)计算),87(),85(),83(),81(f f f f 从而得到:13899.3)]87()85()83()81([812148=++++=f f f f T T ,,14159.334482=-=T T S,14059.31516242=-=S S C .1458.36364121=-=C C R(5)再计算),1615(),1613(),1611(),169(),167(),165(),163(),161(f f f f f f f f从而得到: 14094.316=T30141598=S ,,14159.3,14159.324==R C 51210||-≤-R R , 所以12043.14159.1dx x ≈+⎰3.3高斯求积公式的应用例7 用两点复化高斯求积公式计算10,x I e dx =⎰要求允许误差.106-=ε解 在本算法中取21=+n 时,,110==A A 其中;,)(mab h e x f x -== =++--=∑=)22(2201j jj b a x a b f A a b G.87189637800.1][21)32121()32121(=++-eem =2时, h =21, ]4121)21([4120202j i j j x i f A G +⨯-=∑∑==.57182571650.1)(41341333413341333413=+++=++--eeee m =3时, h =31. .37182769352.1]631)21([6130203=+⨯-=∑∑==j i j j x i f A G .101027.71||||56323--<⨯≈+-G G G3.4 几种方法的比较分析例8 计算积分211ln 2dx x =⎰,精确到0.001.(1)利用矩形公式计算, 因为对于x x f 1)(=,有320()2f x x''<=<(如果1<x <2),所以按照公式0)2(S =+-dx b a x ba . 0<n R <2112n . 如果取n =10,则我们公式的余项的余数得31010.84101200R -<<⨯,我们还必须加进由于在计算函数值实行四舍五入所产生的误差的界限相差于0.16⨯310-,为了这个目的只要计算1x的值到四位小数精确到0.00005就够了.我们有1232527292132152172192 1.051.151.251.351.551.651.751.851.95x x x x x x x x x =========5128.05405.05714.06061.06897.07407.08.08696.09524.02192172152132927252321=========y y y y y y y y y 和6.928469284.0109284.6=(2) 按照梯形公式作同样的计算,在这种情况下,作公式 210,||6n n R R n<<在这儿也试一试取n =10,虽然此时仅可以证3107.16001||-⨯<<n R ,纵坐标是9.18.17.16.15.14.13.12.11.1987654321=========x x x x x x x x x 5263.05556.05882.06250.06667.07143.07692.08333.09091.0987654321=========y y y y y y y y y和1877.669377.01877.621500101=+)( (3) 用辛普森公式做同样的计算作公式 .0))(()2(180)()4(45<≤≤⨯--=n n R b a f n a b R ξξ 并且n =5时有55104.1||-⨯<R .实行计算到五位数字,精确到0.0000058.16.14.12.14321====x x x x 45636.555556.062500.071429.083333.04321和====y y y y 9.17.15.13.11.12927252321=====x x x x x83820.1352632.058824.066667.076923.090909.029********和=====y y y y y.20.150==x x 50000.150000.060000.150和==y y6931525.083820.345636.550000.1301=++)(. 由此可见,用辛普森公式计算得到的值误差最小,计算量相对一般;而用矩形公式计算得到的值误差较大,计算量也比较大;用梯形公式计算的值误差比用矩形公式得到的值要误差小,计算量也是如此.所以我们计算定积分时用辛普森公式往往得到的值误差小,而对没有要求误差大小的,则可以选择辛普森或者是梯形公式,因为这两种方法计算量相对较小.结 束 语本文只讨论了一些一维数值积分方法及其它们的应用,误差分析等有关内容.其中最常用的方法是插值型积分以及复化方法、龙贝格积分方法和高斯积分方法,并讨论了相关求积方法的代数精度和误差分析,并给出了一些例题,分析各种方法的近似值,得出误差分析最小的近似方法.由于篇幅有限,对于高维数值积分方法本文便不再讨论.参考文献[1] 华东师范大学数学系,数学分析(第一版)[M],北京:高等教育出版社,2001. [2] 李庆阳,关治,白峰杉,数值计算原理(第二版)[M],北京: 清华大学出版社, 2008. [3] 肖筱南,现代数值计算方法(第一版)[M],北京: 北京大学出版社, 1999.[4] 菲赫金格尔茨,微积分学教程(第三版)[M],北京: 高等教育出版社, 2005. [5] 裴礼文,数学分析中的典型问题与方法(第一版)[M] ,北京: 北京大学出版社,2004. [6] 李桂成,计算方法(第三版)[M],北京: 高等教育出版社,2010.[7] Yin Y uezhu ,Yang Zhonglian.Calculating Skillfully the Curve Integral and Surface Integral Type 2 bySymmetry, SCIENCE & TECHNOLOGY INFORMATION ,2008(30)The Approximate Numerical Method of the Definite IntegralAbstract This paper mainly discusses common numerical methods of unary function, such as approximate calculation method of interpolation integral, Lebesgue integral and Gauss integration. With these methods in calculating the integral, it will produce some error. In order to reduce the error, we can use after the formula for product and after the Gauss formula. This paper focus on these methods introducing formula of introduction and truncation errors .In addition they can provide examples to analysis size of the error and computation.Keywords interpolation integral Lebesgue integral Gauss integral error analysis approximate computation。
定积分的定义公式分割近似求和取极限
定积分的定义公式分割近似求和取极限定积分这玩意儿,在数学里那可是个相当重要的角色。
它的定义公式——分割近似求和取极限,听起来好像挺复杂,但咱们慢慢捋捋,其实也没那么可怕。
我记得有一次,我在课堂上讲定积分的时候,有个学生一脸迷茫地看着我,那小眼神仿佛在说:“老师,这都是啥呀?”我就跟他说:“别着急,咱们一步一步来。
”咱先说分割。
这就好比你有一块大蛋糕,你要把它切成好多小块。
比如说,一个函数的区间[a,b] ,咱把它分成 n 个小区间,这就是分割。
每个小区间的长度不一定相等,但加起来就是整个区间的长度。
然后是近似。
这就像你切完蛋糕,要估计每一小块的大小。
对于每个小区间里的函数值,咱找个简单的数来近似代替,比如说区间里某一点的函数值。
再说说求和。
把每个小区间里近似的函数值乘以小区间的长度,然后加起来,这就是求和。
最后是取极限。
当把区间分得越来越细,小区间的数量越来越多,每个小区间的长度越来越小,这个求和的结果就会越来越接近一个确定的值,这个值就是定积分的值。
比如说,你要计算从 0 到 1 区间上 x²的定积分。
咱先把这个区间分成 n 个小区间,每个小区间的长度就是 1/n 。
然后在每个小区间里,咱用区间中点的函数值来近似代替。
比如第 i 个小区间的中点是 i/n ,那这个小区间里的函数值就近似为 (i/n)²。
把每个小区间的近似值乘以小区间长度 1/n 再加起来,得到一个式子。
最后让 n 趋向于无穷大,取这个式子的极限,就能得到定积分的值 1/3 。
在实际生活中,定积分也有很多用处呢。
就像你要计算一个不规则图形的面积,或者计算一个物体在一段时间内移动的路程,都能用到定积分。
还记得有一次我装修房子,要计算一面墙的不规则形状的面积,来确定需要多少壁纸。
我就用定积分的思路,把那面墙的形状分割成好多小部分,近似计算每一部分的面积,最后求和取极限,算出了差不多准确的面积,成功买到了合适数量的壁纸。
实验二 怎样计算Pi
数学实验实验报告学院:数学与统计学院班级:数学与应用数学3班学号:0314姓名:康萍时间:实验二怎样计算一、实验目的分别用下列三种方法计算π的近似值,并比较三种方法的精确度: 数值积分法:通过使用编写梯形公式和辛普森公式的程序语言计算π。
泰勒级数法:利用反正切函数泰勒级数计算π。
蒙特卡罗(Monte Carlo )法:通过使用编写蒙特卡罗公式的程序语言来计算π。
二、实验环境基于Windows 环境下的软件。
三、实验的基本理论和方法1、数值积分法以单位圆的圆心为原点建立直角坐标系,则单位圆在第一象限内的部分G 是一个扇形,由曲线])1,0[(12∈-=x x y 及两条坐标轴围成,它的面积4π=S 。
算出了S 的近似值,它的4倍就是π的近似值。
而扇形面积S 实际上就是定积分4112π=-⎰dx x 。
与π有关的定积分有很多,比如211x +的定积分411102π=+⎰dx x 就比21x -的定积分更容易计算,更适合于用来计算π。
一般地,要计算定积分()dx x f ba ⎰,也就是计算曲线()x f y =与直线b x a x y ===,,0所围成的曲边梯形G 的面积S 。
为此,用一组平行于y 轴的直线()b x x x x x a n i x x n n i =<<<<<=-≤≤=-1210,11 将曲边梯形T 分成n 个小曲边梯形,总面积S 分成这些小曲边梯形的面积之和。
如果取n 很大,使每个小曲边梯形的宽度都很小,可以将它上方的边界()()i i x x x x f ≤≤-1近似的看作直线段,将每个小曲边梯形近似的看作梯形来求面积,就得到梯形公式。
如果更准确些,将每个小曲边梯形的上边界近似的看作抛物线段,就得到辛普森公式。
具体公式如下:梯形公式 设分点11,,-n x x 将积分区间],[b a 分成n 等份,即()n i n a b i a x i ≤≤-+=0,/。
试验定积分的近似计算
实验一特殊函数与图形一、问题背景与实验目的著名的Riemann函数大家都很熟悉了,但是关于它的图像你是否清楚呢?除了最上面那几点,其他都很难画吧?你想不想看看下面那些“挤在一起”的点是怎样分布的呢?还有几何中的马鞍面、单叶双曲面等是怎样由直线生成的,是不是也想目睹一下呢?这些,都离不开绘图.实际上绘图一直是数学中的一种重要手段,借助图形,往往可以化繁为简,使抽象的对象得到明白直观的体现.比如函数的基本性质,一个图形常可以使之一目了然,非常有效.它虽不能代替严格的分析与证明,但在问题的研究过程中,可以帮助研究人员节约相当一部分精力.此外,它还可以使计算、证明、建模等的结果得到更明白易懂的表现,有时,这比科学论证更有说服力.同时,数学的教学与学习过程也离不开绘图.借助直观的图形,常可以使初学者更容易接受新知识.如数学分析中有不少函数,其解析式着实让人望而生畏,即使对其性质作了详尽的分析,还是感到难明就里;但如果能看到它的图形,再配合理论分析,则问题可以迎刃而解.又如在几何的学习中,会遇到大量的曲线与曲面,也离不开图形的配合.传统的手工作图,往往费力耗时,效果也不尽理想.计算机恰恰弥补了这个不足,使你可以方便地指定各种视角、比例、明暗,从各个角度进行观察.本实验通过对函数的图形表示和几个曲面(线)图形的介绍,一方面展示它们的特点,另一方面,也将就Matlab软件的作图功能作一个简单介绍.大家将会看到,Matlab 的作图功能非常强大.二、相关函数(命令)及简介1.平面作图函数:plot,其基本调用形式:plot(x,y,s)以x作为横坐标,y作为纵坐标.s是图形显示属性的设置选项.例如:x=-pi:pi/10:pi;y=sin(x);plot(x,y,'--rh','linewidth',2,'markeredgecolor','b','markerfacecolor','g')图1在使用函数plot时,应当注意到当两个输入量同为向量时,向量x与y必须维数相同,而且必须同是行向量或者同是列向量.绘图时,可以制定标记的颜色和大小,也可以用图形属性制定其他线条特征,这些属性包括:linewidth 指定线条的粗细.markeredgecolor 指定标记的边缘色markerfacecolor 指定标记表面的颜色.markersize 指定标记的大小.若在一个坐标系中画几个函数,则plot的调用格式如下:plot(x1,y1,s1,x2,y2,s2,……)2.空间曲线作图函数:plot3,它与plot相比,只是多了一个维数而已.其调用格式如下:plot3(x,y,z,s).例如:x=0:pi/30:20*pi;y=sin(x);z=cos(x);plot3(x,y,z)得到三维螺旋线:图23.空间曲面作图函数:(1)mesh函数.绘制彩色网格面图形.调用格式:mesh(z),mesh(x,y,z)和mesh(x,y,z,c).其中,mesh(x,y,z,c)画出颜色由c指定的三维网格图.若x、y均为向量,则length(x)=n,length(y)=m,[m,n]=size(z).(2)surf在矩形区域内显示三维带阴影曲面图.调用格式与mesh类似.(3)ezmesh用符号函数作三维曲面网格图.调用格式:ezmesh(x,y,z)其中x = x(s,t), y = y(s,t),z = z(s,t).画图区域默认为:-2*pi < s < 2*pi 且-2*pi < t < 2*pi.或者用格式:ezmesh(x,y,z,[smin,smax,tmin,tmax])(4)ezsurf用符号函数作三维曲面图.调用格式与ezmesh类似.(5)sphere画球体命令.4.meshgrid,调用格式:[x,y]=meshgrid(m,n),这里的m,n为给定的向量,可以定义网格划分区域和划分方法.矩阵x和矩阵y是网格划分后的数据矩阵.5.图像的修饰与其他函数:(1)axis equal 控制各个坐标轴的分度,使其相等;(2)colormap设置绘图颜色.调用格式:colormap([r g b])其中r,g,b都是0-1之间的数.或者用格式:colormap(s)s为颜色映像.下面举几个常用的例子:(3)grid网格函数grid on添加网格.grid off取消网格.(4)find找出符合条件的元素在数组中的位置.调用格式:y=find(条件)例如:输入:a=[4 5 78 121 4 665 225 4 1];b=find(a>7)输出:b =3 4 6 7三、实验内容数学分析中,特别是积分部分,我们接触了不少有趣的函数,由于其中有的不是一一对应的,用上面的方法无法画出它们的图像,这时就只能用参数了.此外还有些图形只能用参数来画,比如空间曲线,在计算机上不接受“两个曲面的交线”这种表示,所以也只能用参数来实现.用参数方式作图的关键在于找出合适的参数表示,尤其是不能有奇点,最好也不要用到开方.所以要找的参数最好是有几何意义的.当然这也不可一概而论,需要多积累经验.1.利用函数plot在一个坐标系中画以下几个函数图像,要求采用不同颜色、不同线形、不同的符号标记.函数为:===∈.x t y t z t tπsin(),cos(),sin(2),(0,2)程序如下:t=0:pi/20:2*pi;x=sin(t);y=cos(t);z=sin(2*t);plot(t, x, '--k*', t, y, '-rs', t, z, ':bo')图像如下:图32.绘制类似田螺线的一条三维螺线(方程自己设计).程序如下:t=0:.1:30;x=2*(cos(t)+t.*sin(t));y=2*(sin(t)-t.*cos(t));z=1.5*t;plot3(x,y,-z) %取–z 主要是为了画图看起来更清楚axis equal图像如下:图43.利用函数z=绘制一个墨西哥帽子的图形.程序如下:[a,b]=meshgrid(-8:.5:8); %先生成一个网格c=sqrt(a.^2+b.^2)+eps;z=sin(c)./c;mesh(a,b,z)axis square图像如下:图5思考:这里的 eps 是什么?其作用是什么?4.利用surf 绘制马鞍面图形(函数为:2294x y z =-). 程序如下:[x,y]=meshgrid(-25:1:25,-25:1:25); z=x.^2/9-y.^2/4; surf(x,y,z) title('马鞍面') grid off图像如下:图65.分别用ezmesh 和ezsurf 各绘制一个圆环面,尝试将两个圆环面放在一个图形界面内,观察它们有什么不同之处.提示:圆环面的方程为: 2 ,6 ,)(22222===+-+r R r z R y x ,而圆环面的参数方程为:]2,0[ ],2,0[ ,sin sin )cos (cos )cos (ππ∈∈⎪⎩⎪⎨⎧=+=+=v u u r z v u r R y v u r R x 程序参见附录1. 图像如下:图76.绘制黎曼函数图形,加深对黎曼函数的理解.说明:黎曼函数的定义为1(0,1) 01[01]p p p q x qq q y x x ⎧=∈⎪=⎨⎪=∈⎩,当、为正整数,为既约分数,0,当,及无理点,, 程序参见附录2. 图像如下:图8四、自己动手1. 作出下图所示的三维图形:图9提示:图形为圆环面和球面的组合.2.作出下图所示的墨西哥帽子及其剪裁图形:图10 3.画出球面、椭球面、双叶双曲面、单叶双曲面.4.若要求田螺线的一条轴截面的曲边是一条抛物线:0y时25==.试重新x z设计田螺线的参数方程,并画出该田螺线.5.作出下图所示的马鞍面(颜色为灰色,并有一个标题:“马鞍面”):图116.绘制图8所示的黎曼函数图形,要求分母的最大值n的数值由键盘输入(提示:使用input语句).五、附录附录1:(fulu1.m)程序如下:subplot(1,2,1)ezmesh('(6+2*cos(u))*cos(v)','(6+2*cos(u))*sin(v)','2*sin(u)',[0,2*pi,0,2*pi]) axis equalsubplot(1,2,2)ezsurf('(6+2*cos(u))*cos(v)','(6+2*cos(u))*sin(v)','2*sin(u)',[0,2*pi,0,2*pi])axis equal附录2:(fulu2.m)程序如下:n=100;x=[];y=[];k=1;for i=2:nfor j=1:i-1if gcd(i,j)==1 %用函数gcd(m,n)可求m和n的最大公约数x(k)=j/i;y(k)=1/i;k=k+1;endendendplot(x,y,'.b');axis([0,1,0,1])。
实验二:定积分的近似计算
实验二:定积分的近似计算实验目标:1、 熟悉MATLAB 的程序设计方法;2、 熟悉MATLAB 下命令文件和函数文件的建立和使用;3、 学习定积分的三种近似计算方法:矩形法、梯形法、辛普生法;4、 理解数值计算的误差分析。
问题背景:求定积分的近似值的数值方法就是用被积函数的有限个抽样值的离散或加权平均近似值代替定积分的值。
求某函数的定积分时,在多数情况下,被积函数的原函数很难用初等函数表达出来, 因此能够借助牛顿-莱布尼兹公式计算定积分的机会是不多的。
另外,许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数,对这类函数的定积分,也不能用不定积分方法求解。
由于以上原因,数值积分的理论与方法一直是计算数学研究的基本课题。
定积分的近似解的简单方法包括:矩形公式、梯形公式和辛普生公式。
根据定积分的定义,每个积分和都可以看做是定积分的近似值,即⎰∑=∆=ba ni i i x f dx x f 1)()(ζ在几何意义上,这是用一系列小块区域的面积近似小曲边梯形的面积。
当然,只有当积分区间被分割的很细时,计算结果才具有一定的精确度。
● 矩形法:设积分区间被等分为若干份,第i 份是由][1+i i x x 表示,则该小块区域面积为:)(*1i i i x f x x -+ 或)(*11++-i i i x f x x 或)(*211++-i i i x f x x● 梯形法:设积分区间被等分为若干份,第i 份是由][1+i i x x 表示,取)(i x f 和)(1+i x f 的加权平均值作为平均高度)(i f ζ,则该小块区域面积为:2)()(*11+++-i i i i x f x f x x ● 辛普生法:设积分区间被等分为若干份,第i 份是由][1+i i x x 表示,中点为21+i x ,取函数)(x f 在i x ,1+i x ,21+i x 这是三点的函数值的加权平均值作为平均高度的近似值,则该小块区域面积为:6)()(4)(*1211+++++-i i i i i x f x f x f x x实验内容:1、 试推导定积分的三种近似计算方法的迭代公式(矩形法、梯形法、辛普生法)。
数学分析10.6定积分的近似计算
第十章 定积分的应用 6 定积分的近似计算根据定积分的定义,每一个积分和都可看做是定积分的一个近似值, 如⎰ba f(x )dx=i n 1i i x △)f(x ∑=(或i n1i 1-i x △)f(x ∑=). 这种用一系列小矩形面积来近似表示曲边梯形面积的方法称为矩形法.只有当积分区间被分割得很细时,矩形法才有一定的精确度.一、梯形法将积分区间[a,b]作n 等分,分别依次为a=x 0<x 1<x 2<…<x n <b, △x i =na-b . 相应的被积函数记为y 0,y 1,y 2,…,y n (y i =f(x i ), i=0,1,2,…,n), 并记曲线y=f(x)上相应的点为P 0,P 1,P 2,…,P n (P i (x i ,y i ), i=0,1,2,…,n).将曲线上每一段弧⌒P i-1P i 用弦i 1-i P P 来替代,使得每个小区间[x i-1,x i ]上的曲边梯形换成了真正的梯形,其面积为:2y y i1-i +△x i , i=0,1,2,…,n. 于是,各小梯形面积之和就是曲边梯形面积的近似值,即⎰baf(x )dx=i n1i i 1-i x △2y y ∑=+,亦即⎰b a f(x )dx=n a -b (2y 0+y 1+y 2+…+y n-1+2y n ). 此近似式称为定积分的梯形法公式.二、抛物线法将积分区间[a,b]作n 等分,分别依次为a=x 0<x 1<x 2<…<x 2n <b, △x i =2na-b . 相应的被积函数记为y 0,y 1,y 2,…,y 2n (y i =f(x i ), i=0,1,2,…,2n), 曲线y=f(x)上相应的点为P 0,P 1,P 2,…,P 2n (P i (x i ,y i ), i=0,1,2,…,n).现把区间[x 0,x 2]上的曲线y=f(x)用通过三点P 0(x 0,y 0), P 1(x 1,y 1), P 2(x 2,y 2) 的抛物线p 1(x)= α1x 2+β1x+γ1来近似替代,便有⎰2x x f(x)dx ≈⎰20x x 1(x)p dx=⎰+2x x 1121)γ+x βx (αdx=3α1(x 23-x 03)+2β1(x 22-x 02)+γ1(x 2-x 0) =6x -x 02[(α1x 02+β1x 0+γ1)+(α1x 22+β1x 2+γ1)+α1(x 0+x 2)2+2β1(x 0+x 2)+4γ1] =6x -x 02(y 0+y 2+4y 1)=n6a-b (y 0+4y 1+y 2). 同样的,在[x 2i-2,x 2i ]上,用p i (x)= αi x 2+βi x+γi 来近似替代曲线y=f(x), 可得⎰2i 2-2i x x f(x)dx ≈⎰2i2-2i x x i (x)p dx=n6a-b (y 2i-2+4y 2i-1+y 2i ). 按i=1,2,…,n 把这些近似式相加,得:⎰ba f(x )dx=∑⎰=n1i x x 2i2-2i f(x )dx ≈∑=++n1i 2i 1-2i 2-2i )y y 4y (n 6a -b ,即 ⎰baf(x )dx ≈n6a-b [y 0+y 2n +4(y 1+y 3+…+y 2n-1)+2(y 2+y 4+…+y 2n-2)]. 这就是抛物线法公式,也称为辛普森公式.例:分别用三种求定积分近似值的方法求⎰+12x1dx,并与准确值比较. 解:将区间[0,1]十等分,各分点上被积函数的值如下表(取七位小数):1)用矩形法公式计算,得:⎰+102x 1dx ≈101(y 0+y 1+…+y 9)=0.8099;或⎰+12x 1dx ≈101(y 1+y 2+…+y 10)=0.7600. 2)用梯形法公式计算,得:⎰+102x 1dx ≈101(2y 0+y 1+y 2+…+y 9+2y 10)=0.7850. 3)用抛物线法公式计算,得⎰+102x 1dx ≈n 6a-b [y 0+y 10+4(y 1+y 3+…+y 9)+2(y 2+y 4+…+y 8)]=0.7853982.4)通过牛顿—莱布尼茨公式求原函数得准确值为:⎰+102x 1dx =arctan1=4π=0.78539816… 可见,抛物线法得到的结果最接近准确值.习题1、分别用梯形法和抛物线法近似计算⎰21xdx(将积分区间十等分). 解:将区间[1,2]十等分,各分点上被积函数的值如下表(取七位小数):1)用梯形法公式计算,得:⎰21x dx ≈101(2y0+y 1+y 2+…+y 9+2y 10)=0.69377.2)用抛物线法公式计算,得⎰21x dx ≈n6a-b [y 0+y 10+4(y 1+y 3+…+y 9)+2(y 2+y 4+…+y 10)]=0.693150.注:通过牛顿—莱布尼茨公式求原函数得准确值为:⎰21xdx=ln2=0.693147…2、用抛物线法近似计算⎰π0xsinxdx(分别将积分区间二等分、四等分、六等分). 解:当n=2时,⎰π0x sinx dx ≈12π[1+4(π22+3π22)+2·π2]≈1.852211;当n=4时,⎰π0x sinx dx ≈24π[1+4(π8sin 8π+3π8sin 83π+5π8sin 85π+7π8sin 87π) +2(π22+π2+3π22)]≈1.851937; 当n=6时,⎰π0x sinx dx ≈36π[1+4(π12sin 12π+π22+5π12sin 125π+7π12sin127π +3π22+11π12sin 1211π)+2(π3+2π33+π2+4π33+5π3)]≈1.851940; 注:xsinx的原函数不是初等函数,所以不能直接通过牛顿—莱布尼茨公式求定积分.3、下图为河道某一截面图. 试由测得数据用抛物线法求截面面积. 解:河道的截面面积为: S ≈308[4(0.50+1.30+2.00+1.20+0.55)+ 2(0.85+1.65+1.7 5+0.85)]=8.64(m 2).4、下表所列为夏季某一天每隔两小时测得的气温:(1)按积分平均⎰baf(t)a -b 1dt 求这一天的平均气温,其中定积分值由三种近似法分别计算;(2)若按算术平均∑=121i 1-i C 121或∑=121i i C 121求得平均气温,那么它们与矩形法积分平均和梯形法积分平均各有什么联系?简述理由. 解:(1)设平均气温为T ,n=12, a=0, b=24. 矩形法:T ≈121(C 0+C 1+…+C 11)≈28.71或T ≈121(C 1+C 2+…+C 12)≈28.64. 梯形法:T ≈121(2C 0+C 1+C 2+…+C 11+2C 12)=28.675.抛物线法:T ≈361[C 0+C 12+4(C 1+C 3+…+C 11)+2(C 2+C 4+…+C 10)]≈28.67. (2)∵∑=121i 1-i C 121=121(C 0+C 1+…+C 11);∑=121i i C 121=121(C 1+C 2+…+C 12).∴运用矩形法积分平均的两种形式分别与它们相对应,即为近似值.而梯法积分平均则以C 0和C 12的平均值代替∑=121i 1-i C 121中的C 0或∑=121i iC 121中的C 12,所以也是它们的近似值.。
第七讲-定积分的近似计算
quad 举例
例:用 quad 计算定积分:
dx 0 1 x 2
1
解:
>> quad('1./(1+x.^2)',0,1)
>> quad('1./(1+x.^2)',0,1,10e-10) 函数表达式一定要用 单引号 括起来! 涉及的运算一定要用 数组运算!
dblquad
i 1
n
通常我们取
x1 x2 xn
h ba n
点 i [ xi 1, xi ] 可以任意选取,常见的取法有: 左端点 xi 1 ,右端点
xi 和中点 ( xi 1 xi ) / 2 。
中点法
左点法
右点法
左点法、右点法和中点法
步长
xi h (b a) / n xi a ih, i 1,2, n
抛物线法
设过以上三点的抛物线方程为: y = x2 + x + = p1(x)
则在区间 [x0, x2] 上,有
x2
x0
f ( x)dx p1 ( x)dx x ( x2 x )dx
x2
0
x2
x0
x x x
3 2
x2
x2 x0 (y0 4y1 y2 ) 6 ba (y0 4 y1 y2 ) 6n
i 1 n
x2 i 2
f ( x )dx
ba ( y2i 2 4 y2i 1 y2i ) i 1 6n
抛物线法
整理后可得:
baຫໍສະໝຸດ b a f ( x)dx [ y0 y2n 4( y1 y3 y2n1 ) 6n 2( y2 y4 y2n2 )]
Matlab数值实验定积分的近似计算教程
Matlab数值实验定积分的近似计算教程一、问题背景与实验目的二、相关函数(命令)及简介三、实验内容1.矩形法2.梯形法3.抛物线法4.直接应用Matlab命令计算结果四、自己动手一、问题背景与实验目的利用牛顿—莱布尼兹公式虽然可以精确地计算定积分的值,但它仅适用于被积函数的原函数能用初等函数表达出来的情形.如果这点办不到或者不容易办到,这就有必要考虑近似计算的方法.在定积分的很多应用问题中,被积函数甚至没有解析表达式,可能只是一条实验记录曲线,或者是一组离散的采样值,这时只能应用近似方法去计算相应的定积分.本实验将主要研究定积分的三种近似计算算法:矩形法、梯形法、抛物线法.对于定积分的近似数值计算,Matlab有专门函数可用.二、相关函数(命令)及简介1.sum(a):求数组a的和.2.format long:长格式,即屏幕显示15位有效数字.(注:由于本实验要比较近似解法和精确求解间的误差,需要更高的精度).3.double():若输入的是字符则转化为相应的ASCII码;若输入的是整型数值则转化为相应的实型数值.4.quad():抛物线法求数值积分.格式: quad(fun,a,b) ,注意此处的fun是函数,并且为数值形式的,所以使用*、/、^等运算时要在其前加上小数点,即 .*、./、.^等.例:Q = quad('1./(x.^3-2*x-5)',0,2);5.trapz():梯形法求数值积分.格式:trapz(x,y)其中x为带有步长的积分区间;y为数值形式的运算(相当于上面介绍的函数fun)例:计算x=0:pi/100:pi;y=sin(x);trapz(x,y)6.dblquad():抛物线法求二重数值积分.格式:dblquad(fun,xmin,xmax,ymin,ymax),fun可以用inline 定义,也可以通过某个函数文件的句柄传递.例1:Q1 = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)顺便计算下面的Q2,通过计算,比较Q1 与Q2结果(或加上手工验算),找出积分变量x、y的上下限的函数代入方法.Q2 = dblquad(inline('y*sin(x)'), 0, pi, pi, 2*pi) 例2:Q3 = dblquad(@integrnd, pi, 2*pi, 0, pi)这时必须存在一个函数文件integrnd.m:function z = integrnd(x, y)z = y*sin(x);7.fprintf(文件地址,格式,写入的变量):把数据写入指定文件.例:x = 0:.1:1;y = [x; exp(x)];fid = fopen('exp.txt','w'); %打开文件fprintf(fid,'%6.2f %12.8f\n',y); %写入fclose(fid) %关闭文件8.syms 变量1 变量2 …:定义变量为符号.9.sym('表达式'):将表达式定义为符号.解释:Matlab中的符号运算事实上是借用了Maple的软件包,所以当在Matlab中要对符号进行运算时,必须先把要用到的变量定义为符号.10.int(f,v,a,b):求f关于v积分,积分区间由a到b.11.subs(f,'x',a):将 a 的值赋给符号表达式 f 中的 x,并计算出值.若简单地使用subs(f),则将f的所有符号变量用可能的数值代入,并计算出值.三、实验内容1.矩形法根据定积分的定义,每一个积分和都可以看作是定积分的一个近似值,即在几何意义上,这是用一系列小矩形面积近似小曲边梯形的结果,所以把这个近似计算方法称为矩形法.不过,只有当积分区间被分割得很细时,矩形法才有一定的精确度.针对不同的取法,计算结果会有不同,我们以为例(取),(1)左点法:对等分区间,在区间上取左端点,即取,0.78789399673078,理论值,此时计算的相对误差(2)右点法:同(1)中划分区间,在区间上取右端点,即取,0.78289399673078,理论值,此时计算的相对误差(3)中点法:同(1)中划分区间,在区间上取中点,即取,0.78540024673078,理论值,此时计算的相对误差如果在分割的每个小区间上采用一次或二次多项式来近似代替被积函数,那么可以期望得到比矩形法效果好得多的近似计算公式.下面介绍的梯形法和抛物线法就是这一指导思想的产物.2.梯形法等分区间,相应函数值为().曲线上相应的点为()将曲线的每一段弧用过点,的弦(线性函数)来代替,这使得每个上的曲边梯形成为真正的梯形,其面积为,.于是各个小梯形面积之和就是曲边梯形面积的近似值,,即,称此式为梯形公式.仍用的近似计算为例,取,0.78539399673078,理论值,此时计算的相对误差很显然,这个误差要比简单的矩形左点法和右点法的计算误差小得多.3.抛物线法由梯形法求近似值,当为凹曲线时,它就偏小;当为凸曲线时,它就偏大.若每段改用与它凸性相接近的抛物线来近似时,就可减少上述缺点,这就是抛物线法.将积分区间作等分,分点依次为,,对应函数值为(),曲线上相应点为().现把区间上的曲线段用通过三点,,的抛物线来近似代替,然后求函数从到的定积分:由于,代入上式整理后得同样也有……将这个积分相加即得原来所要计算的定积分的近似值:,即这就是抛物线法公式,也称为辛卜生(Simpson)公式.仍用的近似计算为例,取,=0.78539816339745,理论值,此时计算的相对误差4. 直接应用Matlab命令计算结果(1)数值计算方法1:int('1/(1+x^2)','x',0,1) (符号求积分)方法2:quad('1./(1+x.^2)',0,1) (抛物线法求数值积分)方法3:x=0:0.001:1;y=1./(1+x.^2);trapz(x,y) (梯形法求数值积分)(2)数值计算方法1:int(int('x+y^2','y',-1,1),'x',0,2) (符号求积分)方法2:dblquad(inline('x+y^2'),0,2,-1,1) (抛物线法二重数值积分)四、自己动手1.实现实验内容中的例子,即分别采用矩形法、梯形法、抛物线法计算,取,并比较三种方法的精确程度.2.分别用梯形法与抛物线法,计算,取.并尝试直接使用函数trapz()、quad()进行计算求解,比较结果的差异.3.试计算定积分.(注意:可以运用trapz()、quad()或附录程序求解吗?为什么?)4.将的近似计算结果与Matlab中各命令的计算结果相比较,试猜测Matlab中的数值积分命令最可能采用了哪一种近似计算方法?并找出其他例子支持你的观点.5.通过整个实验内容及练习,你能否作出一些理论上的小结,即针对什么类型的函数(具有某种单调特性或凹凸特性),用某种近似计算方法所得结果更接近于实际值?6.学习fulu2sum.m的程序设计方法,尝试用函数 sum 改写附录1和附录3的程序,避免for 循环.。
定积分近似计算误差c++
定积分近似计算误差c++
定积分的近似计算误差是通过泰勒展开式推导得出的。
假设我
们要计算函数f(x)在区间[a, b]上的定积分,我们可以使用泰勒展
开式来近似计算。
泰勒展开式可以表示为:
f(x) = f(a) + f'(a)(x-a) + f''(a)(x-a)^2/2! +
f'''(a)(x-a)^3/3! + ...
我们可以使用泰勒展开式的前几项来近似计算定积分。
假设我
们使用n阶泰勒展开式来近似计算定积分,那么近似值可以表示为:
∫[a, b] f(x)dx ≈ ∑[i=0 to n] f^(i)(a)(b-
a)^(i+1)/(i+1)!
其中f^(i)(a)表示f(x)的i阶导数在点a处的值。
这个近似值
的误差可以通过泰勒展开式的余项来估计。
余项可以表示为:
R_n = (-1)^(n+1) f^(n+1)(ξ)(b-a)^(n+2)/(n+2)!
其中ξ是区间[a, b]上的某个点。
根据泰勒展开式的余项,我
们可以得到定积分的近似计算误差为|R_n|。
在C++中,我们可以编写一个函数来实现定积分的近似计算,并计算误差。
我们可以使用数值积分方法(如梯形法则、辛普森法则等)来进行近似计算,并根据余项公式来计算误差。
具体实现可以参考数值分析和数值计算的相关算法和库函数,如C++标准库中的数值计算库或者一些开源的数值计算库。
总之,定积分的近似计算误差可以通过泰勒展开式的余项来估计,而在C++中可以编写函数来实现定积分的近似计算,并计算误差。
希望这个回答能够帮助到你。
定积分算法
定积分算法一、概述定积分算法是微积分中的重要概念,用于求解曲线与坐标轴所夹的面积。
在数学中,积分是求解无穷小量的和的过程,而定积分则是求解某个区间内函数值的总和。
通过定积分算法,我们可以计算出曲线下面积的精确值,进而应用于各种实际问题的求解。
二、定积分的基本原理2.1 曲线下面积的近似计算定积分可以通过将曲线下面的区域分割成许多矩形,然后计算每个矩形的面积,最后把这些面积进行累加,从而得到曲线下面的面积。
这种方法称为矩形法。
2.2 矩形法的计算公式假设函数f(x)在区间[a,b]上连续,将区间[a,b]分为n个相等的小区间,每个小区间的长度为Δx=(b-a)/n。
在每个小区间内,选取一个任意点xi,将其对应的函数值f(xi)作为该小区间的高,得到一个矩形,其面积为Δx*f(xi)。
将所有矩形的面积进行累加,即可得到曲线下面的面积的近似值。
2.3 无穷小量与极限的引入上述矩形法可以得到曲线下面积的近似值,但对于曲线比较曲折的情况,求解精确值就无能为力了。
为了解决这个问题,需要引入无穷小量和极限的概念。
将每个小区间的长度Δx取得越来越小,使得Δx无限趋近于0,那么我们可以得到无穷小量dx。
类似地,将每个小区间内任意点xi取得越来越接近于曲线上的某个点x,使得xi无限趋近于x,我们可以得到无穷小量Δx无穷小量dx的关系。
2.4 定积分的计算公式根据无穷小量和极限的引入,可以得到定积分的计算公式:∫(a到b) f(x)dx = lim(n->无穷大) [Σ(f(xi) * Δx)]其中,Σ表示求和,xi为小区间内任意点,Δx为小区间的长度。
三、定积分算法的具体实现3.1 矩形法的算法步骤1.输入函数f(x)。
2.输入区间[a,b]。
3.输入将区间[a,b]分割的份数n。
4.计算Δx=(b-a)/n。
5.根据矩形法的计算公式,进行面积的累加计算。
6.输出曲线下面积的近似值。
3.2 数值积分法的算法步骤数值积分法是通过插值和多项式拟合来求解定积分的近似值。
定积分计算实验报告(3篇)
第1篇一、实验目的1. 理解定积分的概念,掌握定积分的计算方法。
2. 熟悉数值积分的方法,提高数值计算能力。
3. 通过实验,验证定积分的计算结果,加深对定积分理论的理解。
二、实验原理定积分是数学分析中的一个基本概念,它表示函数在某一区间上的累积效果。
对于给定的函数f(x),在区间[a, b]上的定积分可以表示为:∫[a, b] f(x) dx其中,dx表示无穷小的区间宽度。
在实际计算中,定积分往往采用数值积分的方法进行近似计算。
三、实验仪器与软件1. 仪器:计算机2. 软件:MATLAB四、实验步骤1. 输入函数表达式:在MATLAB中输入待积分函数的表达式,例如f(x) = x^2。
2. 设置积分区间:设定积分的上下限a和b。
3. 选择数值积分方法:MATLAB提供了多种数值积分方法,如梯形法、辛普森法、高斯法等。
根据需要选择合适的方法。
4. 进行数值积分计算:调用MATLAB的数值积分函数,如quad函数,进行积分计算。
5. 结果分析:观察计算结果,与理论值进行对比,分析误差来源。
五、实验数据及结果1. 函数表达式:f(x) = x^22. 积分区间:[0, 1]3. 数值积分方法:辛普森法4. 计算结果:I ≈ 1.1666666667六、误差分析1. 理论值:∫[0, 1] x^2 dx = [x^3/3] |[0, 1] = 1/32. 误差来源:a. 数值积分方法的误差:由于数值积分方法是一种近似计算方法,其计算结果与真实值存在一定的误差。
b. 计算过程中的舍入误差:在计算过程中,由于计算机的浮点数表示,可能导致舍入误差。
3. 误差分析:计算结果与理论值相差较大,说明数值积分方法的误差较大。
在实际应用中,可以根据需要选择合适的数值积分方法,以减小误差。
七、实验结论1. 通过本次实验,掌握了定积分的计算方法,了解了数值积分的方法及其优缺点。
2. 了解了数值积分方法在计算过程中的误差来源,为实际应用提供了参考。
求定积分近似值的常用方
求定积分近似值的常用方
积分是数学中用于解释面积的概念,也是求解不定积分的基本工具。
它分成定积分和不定积分,定积分是一种求积分的方法,它的结果是一个常数,是一种完美的积分方法,可以用来计算更复杂的积分。
但定积分有时不容易计算出最终结果,这时可以采用不定积分来近似求解定积分。
求定积分近似值的常用方法共有四种:梯形法、抛物线法、Simpson 法和 Gauss 法。
梯形法是一种定积分近似值计算的最简单最常用的方法。
它是把积分积分区间划分为若干小区间,每个小区间内都保持函数的连续性,把区间[a,b] 分为 n 个小区间,求出每个小区间的积分函数值之和,即可得到定积分的近似值。
抛物线法是采用抛物线去近似求取定积分的值,它也是一种比较简单的方法,其基本思想是把计算的积分区间分成多个小的子区间,将每个小子区间内的函数曲线用抛物线近似代替,并计算抛物线的积分,最后求得积分和,即可得到定积分的近似值。
Simpson 法是由 Simpson 提出的,是一种使用多项式对函数进行拟合近似求解定积分的方法,其基本思想是将定积分积分区间划分为多个子区
间,使用抛物线来近似求取子区间的积分,最后求得抛物线积分的和,即可得到定积分的近似值。
最后,Gauss 法是一种基于多项式拟合近似求定积分的方法,其基本思想是将定积分积分区间划分为多个子区间,使用多项式来近似求取子区间的积分,最后求得多项式积分的和,从而计算出定积分的近似值。
Gauss 法在计算精确度方面要远超过抛物线法和多项式法。
总之,求定积分近似值的常用方法有梯形法、抛物线法、Simpson 法和 Gauss 法。
这四种方法各有优劣,在使用这些方法时应根据不同的需要来选择合适的方法。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
实验二定积分的近似计算一、问题背景与实验目的利用牛顿—莱布尼兹公式虽然可以精确地计算定积分的值,但它仅适用于被积函数的原函数能用初等函数表达出来的情形.如果这点办不到或者不容易办到,这就有必要考虑近似计算的方法.在定积分的很多应用问题中,被积函数甚至没有解析表达式,可能只是一条实验记录曲线,或者是一组离散的采样值,这时只能应用近似方法去计算相应的定积分.本实验将主要研究定积分的三种近似计算算法:矩形法、梯形法、抛物线法.对于定积分的近似数值计算,Matlab有专门函数可用.二、相关函数(命令)及简介1.sum(a):求数组a的和.2.format long:长格式,即屏幕显示15位有效数字.(注:由于本实验要比较近似解法和精确求解间的误差,需要更高的精度).3.double():若输入的是字符则转化为相应的ASCII码;若输入的是整型数值则转化为相应的实型数值.4.quad():抛物线法求数值积分.格式:quad(fun,a,b) ,注意此处的fun是函数,并且为数值形式的,所以使用*、/、^等运算时要在其前加上小数点,即.*、./、.^等.例:Q = quad('1./(x.^3-2*x-5)',0,2);5.trapz():梯形法求数值积分.格式:trapz(x,y)其中x为带有步长的积分区间;y为数值形式的运算(相当于上面介绍的函数fun)例:计算0sin()dx xπ⎰x=0:pi/100:pi;y=sin(x);trapz(x,y)6.dblquad():抛物线法求二重数值积分.格式:dblquad(fun,xmin,xmax,ymin,ymax),fun可以用inline定义,也可以通过某个函数文件的句柄传递.例1:Q1 = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)顺便计算下面的Q2,通过计算,比较Q1 与Q2结果(或加上手工验算),找出积分变量x、y的上下限的函数代入方法.Q2 = dblquad(inline('y*sin(x)'), 0, pi, pi, 2*pi)例2:Q3 = dblquad(@integrnd, pi, 2*pi, 0, pi)这时必须存在一个函数文件integrnd.m:function z = integrnd(x, y) z = y*sin(x);7.fprintf (文件地址,格式,写入的变量):把数据写入指定文件.例:x = 0:.1:1; y = [x; exp(x)];fid = fopen('exp.txt','w'); %打开文件 fprintf(fid,'%6.2f %12.8f\n',y); %写入 fclose(fid) %关闭文件 8.syms 变量1 变量2 …:定义变量为符号. 9.sym('表达式'):将表达式定义为符号.解释:Matlab 中的符号运算事实上是借用了Maple 的软件包,所以当在Matlab 中要对符号进行运算时,必须先把要用到的变量定义为符号. 10.int(f,v,a,b):求f 关于v 积分,积分区间由a 到b .11.subs(f ,'x',a):将 a 的值赋给符号表达式 f 中的 x ,并计算出值.若简单地使用subs(f),则将f 的所有符号变量用可能的数值代入,并计算出值.三、实验内容1. 矩形法根据定积分的定义,每一个积分和都可以看作是定积分的一个近似值,即1()d ()nbi i ai f x x f x ς==∆∑⎰在几何意义上,这是用一系列小矩形面积近似小曲边梯形的结果,所以把这个近似计算方法称为矩形法.不过,只有当积分区间被分割得很细时,矩形法才有一定的精确度.针对不同i ς的取法,计算结果会有不同,我们以 120d 1xx +⎰为例(取100=n ),(1) 左点法:对等分区间b x i n ab a x x a x n i =<<-+=<<<= 10,在区间],[1i i x x -上取左端点,即取1-=i i x ς,12 01d ()1ni i i xf x x ς==∆≈+∑⎰0.78789399673078, 理论值 12 0d 14x x π=+⎰,此时计算的相对误差0.7878939967307840.0031784ππ-=≈(2)右点法:同(1)中划分区间,在区间],[1i i x x -上取右端点,即取i i x =ς,12 01d ()1ni i i xf x x ς==∆≈+∑⎰0.78289399673078, 理论值 12 0d 14x x π=+⎰,此时计算的相对误差0.7828939967307840.0031884ππ-=≈(3)中点法:同(1)中划分区间,在区间1[,]i i x x -上取中点,即取12i ii x x ς-+=, 12 01d ()1ni i i xf x x ς==∆≈+∑⎰0.78540024673078, 理论值 12 0d 14x x π=+⎰,此时计算的相对误差60.7854002467307842.653104ππ--=≈⨯如果在分割的每个小区间上采用一次或二次多项式来近似代替被积函数,那么可以期望得到比矩形法效果好得多的近似计算公式.下面介绍的梯形法和抛物线法就是这一指导思想的产物.2. 梯形法等分区间b x i n a b a x x a x n i =<<-+=<<<= 10,nab x -=∆相应函数值为n y y y ,,,10 (n i x f y i i ,,1,0),( ==). 曲线)(x f y =上相应的点为n P P P ,,,10 (n i y x Pi i i ,,1,0),,( ==) 将曲线的每一段弧i i P P 1-用过点1-i P ,i P 的弦i i P P 1-(线性函数)来代替,这使得每个],[1i i x x -上的曲边梯形成为真正的梯形,其面积为x y y ii ∆⨯+-21,n i ,,2,1 =. 于是各个小梯形面积之和就是曲边梯形面积的近似值,11 11()d ()22nnbi i i i ai i y y x f x x x y y --==+∆≈⨯∆=+∑∑⎰, 即11 ()d ()22bn n ay y b a f x x y y n --≈++++⎰, 称此式为梯形公式.仍用 12 0d 1x x +⎰的近似计算为例,取100=n ,10112 0d ()122n n y y x b a y y x n --≈++++=+⎰ 0.78539399673078, 理论值 12 0d 14x x π=+⎰,此时计算的相对误差60.7853939967307845.305104ππ--=≈⨯很显然,这个误差要比简单的矩形左点法和右点法的计算误差小得多.3. 抛物线法由梯形法求近似值,当)(x f y =为凹曲线时,它就偏小;当)(x f y =为凸曲线时,它就偏大.若每段改用与它凸性相接近的抛物线来近似时,就可减少上述缺点,这就是抛物线法.将积分区间],[b a 作n 2等分,分点依次为b x i n a b a x x a x n i =<<-+=<<<=2102 ,nab x 2-=∆, 对应函数值为n y y y 210,,, (n i x f y i i 2,,1,0),( ==), 曲线上相应点为n P P P 210,,, (n i y x Pi i i 2,,1,0),,( ==). 现把区间],[20x x 上的曲线段)(x f y =用通过三点),(000y x P ,),(111y x P ,),(222y x P 的抛物线)(12x p x x y =++=γβα来近似代替,然后求函数)(1x p 从0x 到2x 的定积分:21 ()d x x p x x =⎰22 ()d x x x x x αβγ++=⎰)()(2)(30220223032x x x x x x -+-+-γβα]4)(2)()()[(62022022202002γβαγβαγβα++++++++++-=x x x x x x x x x x 由于2201x x x +=,代入上式整理后得 21 ()d x x p x x ⎰)](4)()[(612122202002γβαγβαγβα++++++++-=x x x x x x x x )4(621002y y y x x ++-=)4(6210y y y n ab ++-= 同样也有422 ()d x x p x x ⎰)4(6432y y y n ab ++-=……222 ()d n n x nx p x x -⎰)4(621222n n n y y y nab ++-=-- 将这n 个积分相加即得原来所要计算的定积分的近似值:22222212 11()d ()d (4)6ii nnbx i i i i ax i i b af x x p x x y y y n---==-≈=++∑∑⎰⎰,即021******* ()d [4()2()]6bn n n ab af x x y y y y y y y y n---≈++++++++⎰ 这就是抛物线法公式,也称为辛卜生(Simpson )公式.仍用 12 0d 1x x +⎰的近似计算为例,取100=n ,102132124222 0d [4()2()]16n n n x b ay y y y y y y y x n ---≈+++++++++⎰=0.78539816339745,理论值 12 0d 14x x π=+⎰,此时计算的相对误差160.7853981633974542.827104ππ--=≈⨯4. 直接应用Matlab 命令计算结果(1) 数值计算 120d .1xx +⎰ 方法1:int('1/(1+x^2)','x',0,1) (符号求积分)方法2:quad('1./(1+x.^2)',0,1) (抛物线法求数值积分)方法3:x=0:0.001:1; y=1./(1+x.^2);trapz(x,y) (梯形法求数值积分) (2)数值计算 212 01d d x x y y -+⎰⎰方法1:int(int('x+y^2','y',-1,1),'x',0,2) (符号求积分)方法2:dblquad(inline('x+y^2'),0,2,-1,1) (抛物线法二重数值积分)四、自己动手1. 实现实验内容中的例子,即分别采用矩形法、梯形法、抛物线法计算 120d 1xx +⎰,取258=n ,并比较三种方法的精确程度.2. 分别用梯形法与抛物线法,计算 2 1d xx⎰,取120=n .并尝试直接使用函数trapz()、quad()进行计算求解,比较结果的差异.3. 试计算定积分 0sin d xx x+∞⎰.(注意:可以运用trapz()、quad()或附录程序求解吗?为什么?)4. 将 120d 1xx +⎰的近似计算结果与Matlab 中各命令的计算结果相比较,试猜测Matlab 中的数值积分命令最可能采用了哪一种近似计算方法?并找出其他例子支持你的观点.5. 通过整个实验内容及练习,你能否作出一些理论上的小结,即针对什么类型的函数(具有某种单调特性或凹凸特性),用某种近似计算方法所得结果更接近于实际值?6. 学习fulu2sum.m 的程序设计方法,尝试用函数 sum 改写附录1和附录3的程序,避免for 循环.五、附录附录1:矩形法(左点法、右点法、中点法)(fulu1.m ) format long n=100;a=0;b=1;inum1=0;inum2=0;inum3=0; syms x fx fx=1/(1+x^2); for i=1:nxj=a+(i-1)*(b-a)/n; %左点 xi=a+i*(b-a)/n; %右点 fxj=subs(fx,'x',xj); %左点值fxi=subs(fx,'x',xi); %右点值fxij=subs(fx,'x',(xi+xj)/2); %中点值inum1=inum1+fxj*(b-a)/n;inum2=inum2+fxi*(b-a)/n;inum3=inum3+fxij*(b-a)/n;endinum1inum2inum3integrate=int(fx,0,1)integrate=double(integrate)fprintf('The relative error between inum1 and real-value is about: %d\n\n',...abs((inum1-integrate)/integrate))fprintf('The relative error between inum2 and real-value is about: %d\n\n',...abs((inum2-integrate)/integrate))fprintf('The relative error between inum3 and real-value is about: %d\n\n',...abs((inum3-integrate)/integrate))附录2:梯形法(fulu2.m)format longn=100;a=0;b=1;inum=0;syms x fxfx=1/(1+x^2);for i=1:nxj=a+(i-1)*(b-a)/n;xi=a+i*(b-a)/n;fxj=subs(fx,'x',xj);fxi=subs(fx,'x',xi);inum=inum+(fxj+fxi)*(b-a)/(2*n);endinumintegrate=int(fx,0,1)integrate=double(integrate)fprintf('The relative error between inum and real-value is about: %d\n\n',...abs((inum-integrate)/integrate))附录2sum:梯形法(fulu2sum.m),利用求和函数,避免for 循环format longn=100;a=0;b=1;syms x fxfx=1/(1+x^2);i=1:n;xj=a+(i-1)*(b-a)/n; %所有左点的数组xi=a+i*(b-a)/n; %所有右点的数组fxj=subs(fx,'x',xj); %所有左点值fxi=subs(fx,'x',xi); %所有右点值f=(fxi+fxj)/2*(b-a)/n; %梯形面积inum=sum(f) %加和梯形面积求解integrate=int(fx,0,1)integrate=double(integrate)fprintf('The relative error between inum and real-value is about: %d\n\n',...abs((inum-integrate)/integrate))附录3:抛物线法(fulu3.m)format longn=100;a=0;b=1;inum=0;syms x fxfx=1/(1+x^2);for i=1:nxj=a+(i-1)*(b-a)/n; %左点xi=a+i*(b-a)/n; %右点xk=(xi+xj)/2; %中点fxj=subs(fx,'x',xj);fxi=subs(fx,'x',xi);fxk=subs(fx,'x',xk);inum=inum+(fxj+4*fxk+fxi)*(b-a)/(6*n);endinumintegrate=int(fx,0,1)integrate=double(integrate)fprintf('The relative error between inum and real-value is about: %d\n\n',...abs((inum-integrate)/integrate))。