數(shù)值分析方法第六章

1、Chapter 5 數(shù)值微分與數(shù)值積分2先看一個實例:已知20世紀美國人口的統(tǒng)計數(shù)據(jù)為(單位:百萬)年份 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990人口 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4 試計算美國20世紀的年增長率則人口的增長率為時刻的人口為若記),(txt)()(txdtdxtr如何求dtdx 1 數(shù)值微分3 在實際問題中，往往會遇到某函數(shù)在實際問題中，往往會遇到某函數(shù)f(x)f(x)是是用表格用表格表示的表示的, ,用通常的導(dǎo)數(shù)定義無法求導(dǎo)用通常的導(dǎo)數(shù)

2、定義無法求導(dǎo), ,因此要尋求其他因此要尋求其他方法近似求導(dǎo)。常用的數(shù)值微分方法有方法近似求導(dǎo)。常用的數(shù)值微分方法有: :一一. . 運用差商求數(shù)值微分運用差商求數(shù)值微分2 2 運用插值函數(shù)求數(shù)值微分運用插值函數(shù)求數(shù)值微分三三. . 運用樣條插值函數(shù)求數(shù)值微分運用樣條插值函數(shù)求數(shù)值微分四四. . 運用數(shù)值積分求數(shù)值微分運用數(shù)值積分求數(shù)值微分4運用差商求數(shù)值微分運用差商求數(shù)值微分分析分析：微分和積分是一對互逆的數(shù)學(xué)運算。：微分和積分是一對互逆的數(shù)學(xué)運算。復(fù)習(xí)復(fù)習(xí)：導(dǎo)數(shù)：導(dǎo)數(shù) 是差商是差商 當當 時的極限，在精度要求不高的情形下，可時的極限，在精度要求不高的情形下，可以簡單地取差商作為導(dǎo)數(shù)的近似值以

3、簡單地取差商作為導(dǎo)數(shù)的近似值( )fa()( )f ahf ah0h 基本原理：用差商代替微分基本原理：用差商代替微分5可可用用差差商商來來逼逼近近導(dǎo)導(dǎo)數(shù)數(shù)充充分分小小時時當當,h處在點根據(jù)導(dǎo)數(shù)定義x,hhxfhxfhhxfxfhxfhxfxfhhh2)()(lim)()(lim)()(lim)( 0006( )()( )f af ahfah()()( )2f ahf ahfah向后差商近似計算：向后差商近似計算：()( )( )f ahf afah向前差商近似計算：向前差商近似計算：中心差商近似計算：中心差商近似計算： 差商公式差商公式7hxfhxfxf)()()( 000由Taylor展開

4、hxxfhxhfxfhxf002000),( ! 2)( )()(誤差)()( ! 2)()()( )(000hOfhhxfhxfxfxR向前差商8hhxfxfxf)()()( 000由Taylor展開hxxfhxhfxfhxf002000),( ! 2)( )()(誤差)()( ! 2)()()( )(000hOfhhhxfxfxfxR向后差商9hhxfhxfxf2)()()( 000由Taylor展開23000010102300002020()()()()( ),2!3!()()()()(),2!3!hhf xhf xhfxfxfxxhhhf xhf xhfxfxfxhx誤差)()( 6)

5、( )( 12 2)()()( )(22212000hOfhffhhhxfhxfxfxR中心差商10分析分析 1、步長過大則截斷誤差顯著，但如果步長太小又會導(dǎo)、步長過大則截斷誤差顯著，但如果步長太小又會導(dǎo)致舍入誤差的增長。致舍入誤差的增長。 2、在實際計算時，在保證截斷誤差滿足精度要求的前、在實際計算時，在保證截斷誤差滿足精度要求的前提下選取盡可能大的步長提下選取盡可能大的步長 3、事先給出一個合適的步長往往是很困難的。通常在、事先給出一個合適的步長往往是很困難的。通常在在在變步長變步長的過程中實現(xiàn)步長的的過程中實現(xiàn)步長的自動選擇自動選擇。11主要方法：插值多項式，討論函數(shù)主要方法：插值多項式

6、，討論函數(shù)f(x)導(dǎo)數(shù)的近似值求法。導(dǎo)數(shù)的近似值求法。 插值型求導(dǎo)公式插值型求導(dǎo)公式)()(xPxfn 定義：若函數(shù)定義：若函數(shù)f(x)在節(jié)點在節(jié)點 處的函數(shù)處的函數(shù)值已知，作值已知，作f(x)的的n次插值多項式次插值多項式 ，并用，并用 近似代近似代替替f(x)，即，即) 1, 2 , 1(nixi)(xPn由于由于 是多項式，容易求其導(dǎo)數(shù)，故對應(yīng)于是多項式，容易求其導(dǎo)數(shù)，故對應(yīng)于f(x)的每一個插值多項式的每一個插值多項式 ，建立一個數(shù)值微分公式，建立一個數(shù)值微分公式)(xPn)(xPn( )( )nfxP x這樣建立起來的數(shù)值微分公式，稱為這樣建立起來的數(shù)值微分公式，稱為插值型數(shù)值微分公

7、式。)(xPn12(1)(1)11( )( )( )( )( )( )(1)!(1)!nnnnnxfdff xP xxnndx其中其中 之間，上式第一項是之間，上式第一項是 的近的近似值。似值。nxxx,10在( )ifx即使即使 與與f(x)的函數(shù)值處處相差甚微，但兩個的函數(shù)值處處相差甚微，但兩個函數(shù)的導(dǎo)數(shù)在某些點上的值仍可能有很大的差異函數(shù)的導(dǎo)數(shù)在某些點上的值仍可能有很大的差異所以要對誤差進行分析所以要對誤差進行分析)(xPn13125102050100025020015010050300annealing time t, 1000 sSeebeck Coefficient a a, m

8、mV/Ka a 010203040dtda adtda a誤差放大 !)exp()(nKtBAt a a14(1)1( )( )( )( )(1)!nininiffxP xxn如果只計算結(jié)點處的導(dǎo)數(shù)值如果只計算結(jié)點處的導(dǎo)數(shù)值( )( )inifxP x數(shù)值微分公式的余項為數(shù)值微分公式的余項為1510001101( )( )( )( )2, , ( )( )( )( )2f xf xhf xfhx xf xf xhf xfh（n=1）011010110( )()()xxxxp xf xf xxxxx100101( )( )( )( )( )( )f xf xf xhf xf xf xh16200

9、122102220121() 3 ()4 ( )()( )231( )()()( )261() () 4 ( ) 3 ()( )23hf xf xf xf xfhhf xf xf xfhhf xf xf xf xfh02,x x類似地，可以得到類似地，可以得到n=2的帶余項的三點公式的帶余項的三點公式172(4)00121222(4)101222(4)20121221() ()2 ()() ()()61() ()2 ()()( )121() ()2 ()() ()()6hfxf xf xf xhffhhfxf xf xf xfhhfxf xf xf xhffh 02,x x類似地，可以得到類似

10、地，可以得到n=2的帶余項的二階三點公式的帶余項的二階三點公式18利用三次樣條插值函數(shù)構(gòu)造數(shù)值微分公式利用三次樣條插值函數(shù)構(gòu)造數(shù)值微分公式對于給定函數(shù)表對于給定函數(shù)表 和適當?shù)剡吔鐥l件，有三次樣條插值函數(shù)和適當?shù)剡吔鐥l件，有三次樣條插值函數(shù)S(x)，并用并用S(x)近似代替近似代替f(x)，即，即由于由于S(x)是一個分段三次多項式，在各子區(qū)間是一個分段三次多項式，在各子區(qū)間 上容易求出其導(dǎo)數(shù)，故建立一個數(shù)上容易求出其導(dǎo)數(shù)，故建立一個數(shù)值微分公式值微分公式.,)()(baxxSxf,1iixx), 2 , 1(ni192211111( ) ( )()() 2,()6iiiiiiiiiiifxS

11、xM xxMxxhhf xxMM( )iiSxM注意注意1111( ) ( )()()(,1,2,)iiiiiiiifxSxMxxMxxhxxxi數(shù)值微分公式數(shù)值微分公式20例：例： 利用函數(shù)利用函數(shù) 在節(jié)點在節(jié)點 上的函數(shù)值和邊界條件上的函數(shù)值和邊界條件S(-1)=0.0740, S(1)=-0.0740構(gòu)造三次樣條插值函數(shù)構(gòu)造三次樣條插值函數(shù)S(x)，并用它來計算，并用它來計算f(x)和和f(x)在下列點在下列點處的近似值。處的近似值。22511)(xxf1 0.1 (0,1,20)ixi i )100, 1 , 0(02. 01kkxi2122MATLAB中的數(shù)值微分中的數(shù)值微分數(shù)值微分

12、的實現(xiàn)數(shù)值微分的實現(xiàn)在在MATLAB中，沒有直接提供求數(shù)值導(dǎo)數(shù)的函數(shù)，只有中，沒有直接提供求數(shù)值導(dǎo)數(shù)的函數(shù)，只有計算向前差分的函數(shù)計算向前差分的函數(shù)diff，其調(diào)用格式為：，其調(diào)用格式為：DX=diff(X)：計算向量：計算向量X的向前差分，的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,n-1。DX=diff(X,n)：計算：計算X的的n階向前差分。階向前差分。diff(X,2)=diff(diff(X)。23例例 設(shè)設(shè)x由由0,2間均勻分布的間均勻分布的10個點組成，求個點組成，求sinx的的13階差分。階差分。命令如下命令如下：X=linspace(0,2*pi,10);Y=

13、sin(X);DY=diff(Y); %計算計算Y的一階差分的一階差分D2Y=diff(Y,2); %計算計算Y的二階差分的二階差分，也可用命令也可用命令diff(DY)計算計算D3Y=diff(Y,3); %計算計算Y的三階差分的三階差分，也可用也可用diff(D2Y)或或diff(DY,2)24例例 用不同的方法求函數(shù)用不同的方法求函數(shù)f(x)的數(shù)值導(dǎo)數(shù)，并在同一個坐標系中的數(shù)值導(dǎo)數(shù)，并在同一個坐標系中做出做出f(x)的圖像。的圖像。程序如下：程序如下：f=inline(sqrt(x.3+2*x.2-x+12)+(x+5).(1/6)+5*x+2);g=inline(3*x.2+4*x-1

14、)./sqrt(x.3+2*x.2-x+12)/2+1/6./(x+5).(5/6)+5);x=-3:0.01:3;p=polyfit(x,f(x),5); %用用5次多項式次多項式p擬合擬合f(x)dp=polyder(p); %對擬合多項式對擬合多項式p求導(dǎo)數(shù)求導(dǎo)數(shù)dpdpx=polyval(dp,x); %求求dp在假設(shè)點的函數(shù)值在假設(shè)點的函數(shù)值dx=diff(f(x,3.01)/0.01; %直接對直接對f(x)求數(shù)值導(dǎo)數(shù)求數(shù)值導(dǎo)數(shù)gx=g(x); %求函數(shù)求函數(shù)f的導(dǎo)函數(shù)的導(dǎo)函數(shù)g在假設(shè)點的導(dǎo)數(shù)在假設(shè)點的導(dǎo)數(shù)plot(x,dpx,x,dx,.,x,gx,-); %作圖作圖252 數(shù)值

15、積分數(shù)值積分在高等數(shù)學(xué)中，牛頓在高等數(shù)學(xué)中，牛頓-萊布尼茲公式萊布尼茲公式F(x)是被積函數(shù)是被積函數(shù)f(x)的原函數(shù)的原函數(shù)完美地解決了定積分問題。完美地解決了定積分問題。)a(F)b(F)x(Fdx)x(fbaba 26 在工程技術(shù)和科學(xué)研究中，常常遇在工程技術(shù)和科學(xué)研究中，常常遇到如下情況：到如下情況：1. f(x)的結(jié)果復(fù)雜，求原函數(shù)困難。的結(jié)果復(fù)雜，求原函數(shù)困難。2. f(x)的原函數(shù)不存在，或不能用的原函數(shù)不存在，或不能用初等函數(shù)表示。初等函數(shù)表示。3. f(x)無函數(shù)式，只有函數(shù)表。無函數(shù)式，只有函數(shù)表。27R1R2F 顆粒尺寸，顆粒尺寸，R 分布密度，分布密度，f 例：計算 R

16、1 R R2 的顆粒數(shù)量(分數(shù)) 21)(RRdRRfF原始方法：曲線下面面積所占方格數(shù)目數(shù)值積分雛形 數(shù)格子28當當 f (x) 是閉區(qū)間上的連續(xù)函數(shù)時，如何是閉區(qū)間上的連續(xù)函數(shù)時，如何計算下列定積分？計算下列定積分？badxxffI)()(積分中值定理積分中值定理)()()()(fabdxxffIba定積分的計算定積分的計算29上述公式的幾何意義上述公式的幾何意義底為底為 而高為而高為 的矩形面積恰好等于的矩形面積恰好等于所求曲邊梯形的面積所求曲邊梯形的面積 。ba( )f( )I f問題問題：如何確定：如何確定 ？I(f)Oxy30)()()()(afabdxxffIba)()()()(

17、bfabdxxffIba)2()()()(bafabdxxffIba( )( )( )( )()2baf af bI ff x dxba左矩形左矩形右矩形右矩形矩形矩形中矩形中矩形定積分常用的近似公式定積分常用的近似公式31 線性近似：梯形公式ab梯形近似xyy = f(x)( )( )( )2babaf x dxf af b簡單，但是ab梯形近似xyy = f(x)32 拋物線近似：辛普森(Simpson)公式xyy = f(x)ab2ba 拋物線近似 )(24)(6)(bfbafafabdxxfba通常，拋物線近似優(yōu)于梯形近似但是，有特例 ab2ba xyy = f(x)33利用梯形公式，

18、利用梯形公式，Simposn公式公式公式計算公式計算 ，并與精確值進行比較。，并與精確值進行比較。dxxI15 . 04267767. 0) 15 . 0(25 . 015 . 0dxxI解：梯形公式解：梯形公式Simposn公式公式10.50.5( 0.5 4 0.75 1) 0.430934036Ixdx43096441. 03215 . 0315 . 0 xdxxI精確值精確值34插值型求積公式插值型求積公式 利用插值多項式來構(gòu)造求積公式：在積分利用插值多項式來構(gòu)造求積公式：在積分區(qū)間區(qū)間a, b上取一組點上取一組點 用用f(x)的的n次插值多項式次插值多項式 來近似代替被積函數(shù)來近似代

19、替被積函數(shù)f(x) bxxxan 10 nkkkn)x(l )x(f)x(L00( )()( )nbbkkaakf x dxf xlx dx35若記若記 bankjjjkjbakkdxxxxxdx)x(lA0則得數(shù)值求積公式則得數(shù)值求積公式 nkkkba)x(fAdx)x(f0稱為稱為插值型求積公式插值型求積公式 36數(shù)值積分公式的余項數(shù)值積分公式的余項1(1)1101 ( )()( ) ( )( )( )(1)!( , )( )()()()nbnkkaknbbnnaannRff x dxA f xff xL x dxx dxna bxxxxxxx37數(shù)值積分公式的代數(shù)精度數(shù)值積分公式的代數(shù)精

20、度nkkkbaxfAdxxf0)()(如果對于任意不高于如果對于任意不高于m次的多項式都準確地成立，次的多項式都準確地成立，而對于而對于 不能準確地成立，則稱該求積公式不能準確地成立，則稱該求積公式的的代數(shù)精度為代數(shù)精度為m.1mx數(shù)值求積方法是近似方法，為保證精度，希望所提供數(shù)值求積方法是近似方法，為保證精度，希望所提供求積公式對于求積公式對于“盡可能多盡可能多”的函數(shù)是準確的。的函數(shù)是準確的。38幾個常用的求積公式的代數(shù)精度1.梯形公式的代數(shù)精度222( )1( ) ( ( )( )( )2( )11( )()22 ( ( )( )()( )22bbaababbbaaabaf xf x d

21、xdxbabaT ff af bbaf x dxf xxf x dxxdxxbababaT ff af babf x dx當時當時392233322( )11( )()33 ( ( )( )()( )22bbbaaabaf xxf x dxx dxxbababaT ff af babf x dx當時所以梯形公式具有一次的代數(shù)精度402.辛普森公式的代數(shù)精度2222( )1( )( )4()( )62( )( )1( )()2( )4()( )621(4)()622( )bbaababbaabafxfx dxdxbababaS ff aff bbafx dxS ffxxfx dxxdxbabab

22、aS ff aff bbaababbafx dxS f當時所 以成 立當時所 以成 立4122332222233( )1( )()3 ( ( )4 ()( )62(4()621(222)()63( ) bbaabaf xxf x dxx dxbababaS ff aff bbaababbaaabbbaf x dxS f當時即精確成立423344333332233322344( )1( )()4 ( ( )4 ()( )62(4()621(33)6231()()624( ) bbaabaf xxf x dxx dxbababaS ff aff bbaababbaaaa babbbbaaa bab

23、bbaf x dxS f當時即精確成立可以證明對可以證明對 等式不成立等式不成立辛普森公式具有三次代數(shù)精度辛普森公式具有三次代數(shù)精度4( )f xx43n +1個結(jié)點的插值型數(shù)值積分公式至少有個結(jié)點的插值型數(shù)值積分公式至少有 n 階代數(shù)精度。階代數(shù)精度。(1)( )( )(1)!nbafR fx dxn因此對于次數(shù)不大于因此對于次數(shù)不大于n的多項式的多項式其余項其余項 0R f 44 求插值型求積公式求插值型求積公式 )21()21()(1011fAfAdxxf 并確定其代數(shù)精度。并確定其代數(shù)精度。分析：分析：實際上該題目是求實際上該題目是求A0，A1，并確定其，并確定其代數(shù)精度代數(shù)精度。21

24、,2110 xx因因為為是是插插值值型型的的，且且， 11)21()21()(ffdxxf從而求積公式為從而求積公式為有有2個結(jié)點個結(jié)點,因而代數(shù)精度大于等于因而代數(shù)精度大于等于1，從，從m=2開始驗證代數(shù)精度開始驗證代數(shù)精度。，從從而而1 m解解 1100)(dxxlA 1111)(dxxlA 11232dxx21)21()21( ff,)(2xxf 對對 11, 1)21(dxx1)21(11 dxx45在等分情況下，稱為牛頓在等分情況下，稱為牛頓-科茨公式科茨公式將將a, bn等分等分 為步長，節(jié)點為等分點為步長，節(jié)點為等分點nabh (0,1,2,., )kxakhknx=a+thx

25、xj =(t - j)h4600()()nnkjj ktj hAhdtkj h bankjjjkjbakkdxxxxxdx)x(lA0 x xj =(t - j)hdx =hdt00( )00( 1)()!()!( 1)()()()!()!n knnjj kn knnnkjj khtj dtk nkbatj dtba Cn k nk科茨系數(shù)科茨系數(shù)47 科特斯系數(shù)表48 且滿足 注：注：Cotes 系數(shù)系數(shù)僅取決于僅取決于 n 和和 i，可通過查表得到。與，可通過查表得到。與被積函數(shù)被積函數(shù) f (x) 及積分區(qū)間及積分區(qū)間 a, b 均無關(guān)。均無關(guān)。q 牛頓牛頓- -科特斯公式科特斯公式ni

26、inibaxfCabxxf0)()()(d)(2) )()(ninniCC(1) 10)(niniC49梯形公式梯形公式 當當n=1, x0=a, x1=b時，有時，有 2201011010abdxabaxdxxxxxAabdxbabxdxxxxxAbabababa2)b(f)a(fabdx)x(fba 50辛普森公式辛普森公式 當當n=2, x0=a, x1=(a+b)/2, x2=b時，有時，有6646120210221012012010210abdx)xx)(xx()xx)(xx(Aabdx)xx)(xx()xx)(xx(Aabdx)xx)(xx()xx)(xx(Abababa)b(f)

27、ba(f)a(fabdx)x(fba24651柯特斯公式柯特斯公式其中其中)4 , 3 , 2 , 1 , 0(4kabkaxk01234()( )7 ( ) 32 ( ) 12 ( ) 32 ( ) 7 ( )90bab af x dxf xf xf xf xf x52 定理定理：若：若 在在a,b上連續(xù)，梯形公式的截斷上連續(xù)，梯形公式的截斷誤差為誤差為)( xf)(12)(31fabfR若若 在在a,b上連續(xù)，辛普生公式的截斷誤差為上連續(xù)，辛普生公式的截斷誤差為)()4(xf)()2(901)4(52fabfR若若 在在a,b上連續(xù)，柯特斯公式的截斷誤差為上連續(xù)，柯特斯公式的截斷誤差為)(

28、) 6 (xf7(6)48()( )9454baRff 53實際計算中，當積分區(qū)間較大時，直接使用這些積分公式精度難以保證。為了提高計算精度，采用復(fù)合求積的方法。將區(qū)間a,b適當分割成若干個子區(qū)間，對每個子區(qū)間使用求積公式，構(gòu)成所謂的復(fù)化求積公式，這是提高積分精度的一個常用的方法。54xyy = f(x)ab將 a, b 分成若干小區(qū)間xk-1, xk, (k = 1, 2, , n)x0 x1x2xnxn-1h前提：節(jié)點函數(shù)值 f(xk) 已知11( )( )2()( )2bnkkahf x dxf af xf b復(fù)合梯形公式復(fù)合梯形公式55將積分區(qū)間將積分區(qū)間 a,b 劃分為劃分為n等分等

29、分, ,記子區(qū)間記子區(qū)間 的中點為的中點為 ，在每個小區(qū)間上應(yīng)用辛普森，在每個小區(qū)間上應(yīng)用辛普森公式，則有公式，則有 1,kkxxhxxkk2121)()(4)(6)()(11010211kkkbankxxnkxfxfxfhdxxfdxxfIkk)()(2)(4)(6101121bfxfxfafhnknkkk)()(2)(4)(6101121bfxfxfafhSnknkkkn記記 復(fù)合辛普森公式復(fù)合辛普森公式5611110042113014( )7 ( )32() 12()9032() 14()7 ( )nnbnakkkknnkkkkhf x dxCf af xf xf xf xf b復(fù)合柯特

30、斯公式復(fù)合柯特斯公式hxxhxxhxxkkkkkk43;21;41432141其中其中57誤差分析誤差分析A、復(fù)合梯形公式的誤差、復(fù)合梯形公式的誤差)(12)( )( 12)(22fhabafbfhIdxxfnba或B、復(fù)合、復(fù)合Simpson公式的誤差公式的誤差4(3)(3)4(4 )1( )() ( )( )1802()()1802bnahfx dxSfbfabahf 或58C、復(fù)合、復(fù)合Cotes公式的誤差公式的誤差6(5 )(5 )6( 6 )2()() ( )()94542()()()9454bnahfx dxCfbfabahf 或59 用復(fù)化梯形公式計算定積分用復(fù)化梯形公式計算定積

31、分 才能使誤差不超過才能使誤差不超過 10dxeIx51021解解: :取取 , ,則則 , ,又區(qū)間長度又區(qū)間長度b-a=1b-a=1，對，對復(fù)化梯形公式有余項復(fù)化梯形公式有余項 xexf)(xexf )(52210211121)(12)( enfhabxRT即即 ,n212.85,n212.85，取，取n=213n=213，即將區(qū)間，即將區(qū)間0,10,1分為分為213213等份時，用復(fù)化梯形公式計算誤差等份時，用復(fù)化梯形公式計算誤差不超過不超過 。 52106en51021問區(qū)間問區(qū)間0,10,1應(yīng)分多少等份應(yīng)分多少等份60p復(fù)合求積方法對提高精度是一種行之有效的方法，但是復(fù)合求積方法對提

32、高精度是一種行之有效的方法，但是在使用求積公式之前必須先給出步長。在使用求積公式之前必須先給出步長。p步長取得太大精度難以保證，步長太小則會導(dǎo)致計算量步長取得太大精度難以保證，步長太小則會導(dǎo)致計算量的增加，而事先給出一個合適的步長往往是困難的。的增加，而事先給出一個合適的步長往往是困難的。p為了解決這個問題，較好的方法就是采用為了解決這個問題，較好的方法就是采用變步長變步長的計算的計算方法方法p即在步長逐次減半的過程中，反復(fù)利用復(fù)合求積公式進即在步長逐次減半的過程中，反復(fù)利用復(fù)合求積公式進行計算，直到所求得的積分值滿足精度要求為止。行計算，直到所求得的積分值滿足精度要求為止。610123456

33、7891000.511.501234567891000.511.501234567891000.511.562 設(shè)將積分區(qū)間設(shè)將積分區(qū)間 a,b n等分，即分成等分，即分成n個子區(qū)間，個子區(qū)間，一共有一共有n+1個節(jié)點，即個節(jié)點，即x=a+kh, , k=0,1,，n，步，步長長 。對于某個子區(qū)間。對于某個子區(qū)間 , ,利用梯形公利用梯形公式計算積分近似值有式計算積分近似值有 nabh1,kkxx)()(21kkxfxfh11011()()2( )2()( )2nnkkknkkhTf xf xhf af xf b對整個區(qū)間對整個區(qū)間a,ba,b有有63將子區(qū)間將子區(qū)間 再二等份再二等份, ,取

34、其中點取其中點作新節(jié)點作新節(jié)點, ,此時區(qū)間數(shù)增加了一倍為此時區(qū)間數(shù)增加了一倍為2n,2n,對某個子區(qū)對某個子區(qū)間間 , ,利用復(fù)化梯形公式計算其積分近似值利用復(fù)化梯形公式計算其積分近似值 。1,kkxx)(21121kkkxxx1,kkxx)()(2)(4121kkkxfxfxfh對整個區(qū)間對整個區(qū)間a,ba,b有有 1012)()(2)(421nkkkknxfxfxfhT10101)(2)()(421nkknkkkxfhxfxfh比較比較 和和 有有nTnT2102)(2221nkknnxfhTT64 當把積分區(qū)間分成當把積分區(qū)間分成n等份，用復(fù)化梯形等份，用復(fù)化梯形公式計算積分公式計算積

35、分I的近似值的近似值 時，截斷誤差為時，截斷誤差為 nT)(122nnnfnababTIR 若把區(qū)間再分半為若把區(qū)間再分半為2n2n等份，計算出定積分等份，計算出定積分的近似值的近似值 ，則截斷誤差為，則截斷誤差為 nT2)(2122222nnnfnababTIR 當當 在區(qū)間在區(qū)間a,ba,b上變化不大時上變化不大時, ,有有 )(xf )()(2nnff 214nnITIT所以所以 65 可見可見, ,當步長二分后誤差將減至當步長二分后誤差將減至 , ,將將上式移項整理，可得驗后誤差估計式上式移項整理，可得驗后誤差估計式 41)(3122nnnTTTI上式說明，只要二等份前后兩個積分值上式

36、說明，只要二等份前后兩個積分值和和 相當接近，就可以保證計算結(jié)果相當接近，就可以保證計算結(jié)果的誤差很小，使的誤差很小，使 接近于積分值接近于積分值I I。 nTnT2nT2nT266變步長的梯形求積算法實現(xiàn)變步長的梯形求積算法實現(xiàn)（1 1）變步長的梯形求積法的計算步驟）變步長的梯形求積法的計算步驟 變步長梯形求積法。它是以梯形求積公式為變步長梯形求積法。它是以梯形求積公式為基礎(chǔ)，逐步減少步長，基礎(chǔ)，逐步減少步長，按如下遞推公式求二分后的按如下遞推公式求二分后的梯形值梯形值102)(2221nkknnxfhTT其中其中Tn和和T2n分別代表二等分前后的積分值分別代表二等分前后的積分值 如果如果

37、, (為給定的誤差限為給定的誤差限 ) 則則T2n作為積分的近似值作為積分的近似值, 否則繼續(xù)進行二等分否則繼續(xù)進行二等分, 即即轉(zhuǎn)轉(zhuǎn) 再計算，直到滿足所要求的精度為止，最終取再計算，直到滿足所要求的精度為止，最終取二分后的積分值二分后的積分值T2n 作為所求的結(jié)果作為所求的結(jié)果 nnTT2nnTThh2,267例例 用變步長梯形求積法計算定積分用變步長梯形求積法計算定積分解解: : 先對整個區(qū)間先對整個區(qū)間 0,10,1 用梯形公式用梯形公式, ,對于對于 10dsinxxxI8410709. 0) 1 (, 1)0(,sin)(ffxxxf所以有所以有 9207355. 0) 1 ()0(

38、211ffT然后將區(qū)間二等份然后將區(qū)間二等份, ,由于由于 , ,故有故有 12( )0.958851f9397933.0)(21212112fTT進一步二分求積區(qū)間進一步二分求積區(qū)間, ,并計算新分點上的函數(shù)值并計算新分點上的函數(shù)值 9088516.0)(,9896158.0)(4341ff68有有 9445135.0)()(4121434124ffTT這樣不斷二分下去，計算結(jié)果如這樣不斷二分下去，計算結(jié)果如P P117117所示。所示。積分的準確值為積分的準確值為0.9460830.946083。 69龍貝格求積公式龍貝格求積公式 變步長梯形求積法算法簡單，但精度較差，收斂變步長梯形求積法

39、算法簡單，但精度較差，收斂速度較慢，但可以利用梯形法算法簡單的優(yōu)點，形速度較慢，但可以利用梯形法算法簡單的優(yōu)點，形成一個新算法，這就是龍貝格求積公式。龍貝格公成一個新算法，這就是龍貝格求積公式。龍貝格公式又稱逐次分半加速法。式又稱逐次分半加速法。 根據(jù)積分區(qū)間分成根據(jù)積分區(qū)間分成n等份和等份和2n等份時的誤差估計等份時的誤差估計式可得式可得 )(3122nnnTTTI 積分值積分值 的誤差大致等于的誤差大致等于 , ,如果用如果用 對對 進行修正時，進行修正時， 與與 之和比之和比 更接近積分真值更接近積分真值nT2)(312nnTT)(312nnTTnT2)(312nnTTnT2nT2702

40、22141()333nnnnnnTTTTTT考察考察 與與n等份辛普森公式等份辛普森公式 之間的關(guān)系。之間的關(guān)系。將復(fù)化梯形公式將復(fù)化梯形公式 nS)()(2)(211bfxfafhTnkkn梯形變步長公式梯形變步長公式 102)(2221nkknnxfhTT代入代入 表達式得表達式得 nT121100( )4()2()( )6nnnknkkkhTf af xf xf bSnnnTTS31342故故 nT71nnnnnTTTTTT3134)(31222 (6.9) TnS這就是說，用梯形法二分前后兩個積分值這就是說，用梯形法二分前后兩個積分值 和和 作線性組合，結(jié)果卻得到復(fù)化辛普森公式計算得到

41、作線性組合，結(jié)果卻得到復(fù)化辛普森公式計算得到的積分值的積分值 。 nTnT2nS72再考察辛普森法。其截斷誤差與再考察辛普森法。其截斷誤差與 成正比，因此，成正比，因此，如果將步長折半，則誤差減至如果將步長折半，則誤差減至 ，即有，即有 4h1/161612nnSISI由此可得由此可得 nnSSI15115162可以驗證可以驗證, ,上式右端的值其實等于上式右端的值其實等于C Cn n，就是說，用辛，就是說，用辛普森公式二等份前后的兩個積分值普森公式二等份前后的兩個積分值S Sn n和和S S2n2n 作線性組作線性組合后，可得到柯特斯公式求得的積分值合后，可得到柯特斯公式求得的積分值C Cn

42、 n，即有，即有 nnnSSC1511516273用同樣的方法，根據(jù)柯特斯公式的誤差公式，可進用同樣的方法，根據(jù)柯特斯公式的誤差公式，可進一步導(dǎo)出龍貝格公式一步導(dǎo)出龍貝格公式 nnnCCR63163642p在變步長的過程中運用上述公式，就能將粗糙的梯在變步長的過程中運用上述公式，就能將粗糙的梯形值形值Tn逐步加工成精度較高的辛普森值逐步加工成精度較高的辛普森值Sn、柯特斯值、柯特斯值Cn和龍貝格值和龍貝格值Rnp或者說，將收斂緩慢的梯形值序列或者說，將收斂緩慢的梯形值序列Tn加工成收斂迅加工成收斂迅速的龍貝格值序列速的龍貝格值序列Rn，這種加速方法稱為龍貝格算法，這種加速方法稱為龍貝格算法74

43、龍貝格求積法計算步驟龍貝格求積法計算步驟 用梯形公式計算積分近似值用梯形公式計算積分近似值 按變步長梯形公式計算積分近似值按變步長梯形公式計算積分近似值 將區(qū)間逐次分半將區(qū)間逐次分半, ,令區(qū)間長度令區(qū)間長度 )()(21bfafabT), 2 , 1 , 0(2kabhk102)(2221nkknnxfhTT計算計算)2(kn 按加速公式求加速值按加速公式求加速值 322nnnnTTTS1522nnnnSSSC6322nnnnCCCR梯形加速公式：梯形加速公式： 辛普森加速公式：辛普森加速公式： 龍貝格求積公式：龍貝格求積公式： 75 精度控制；直到相鄰兩次積分值精度控制；直到相鄰兩次積分值

44、 nnRR2（其中（其中為允許的誤差限）則終止計算并取為允許的誤差限）則終止計算并取R Rn n作為積分作為積分 的近似值，否則將區(qū)間再對分，重的近似值，否則將區(qū)間再對分，重復(fù)復(fù) ， 的計算，直到滿足精度要求為止。的計算，直到滿足精度要求為止。 badxxf)(76例例 用龍貝格算法計算定積分用龍貝格算法計算定積分 要求相鄰兩次龍貝格值的偏差不超過要求相鄰兩次龍貝格值的偏差不超過解解: :由題意由題意 102d14xxI510214)(, 1,0 xxfba3)24(21) 1 ()0(211ffT1 . 351621321)(21212112fTT13118. 3)56. 2764. 3(4

45、11 . 321)()(4121434124ffTT13899. 3)()()()(81218785838148ffffTT14094. 3)()()()()()()()(16121161516131611169167165163161816ffffffffTT771333.33134121TTS14157. 33134242TTS14159.33134484TTS14159.331348168TTS14212. 31511516121SSC14159. 31511516242SSC14159. 31511516484SSC7814158. 36316364121CCR14159. 36316

46、364242CCR由于由于 ，于是有，于是有 00001. 012RR14159. 3d14102xxI79本章小結(jié)本章小結(jié) 本章介紹了積分和微分的數(shù)值計算方法，其基本本章介紹了積分和微分的數(shù)值計算方法，其基本原理主要是逼近論，即設(shè)法構(gòu)造某個簡單函數(shù)原理主要是逼近論，即設(shè)法構(gòu)造某個簡單函數(shù)P(x)P(x)近近似表示似表示f(x)f(x)，然后對，然后對P(x)P(x)求積或求導(dǎo)得到求積或求導(dǎo)得到f(x)f(x)的積分的積分或?qū)?shù)的近似值。基于插值原理，推導(dǎo)了數(shù)值積分和或?qū)?shù)的近似值。基于插值原理，推導(dǎo)了數(shù)值積分和數(shù)值微分的基本公式。數(shù)值微分的基本公式。80MATLAB數(shù)值積分數(shù)值積分基于變步長

47、辛普生法，基于變步長辛普生法，MATLAB給出了給出了quad函數(shù)來求定積分。函數(shù)來求定積分。該函數(shù)的調(diào)用格式為：該函數(shù)的調(diào)用格式為： I,n=quad(fname,a,b,tol,trace)fname是被積函數(shù)名。是被積函數(shù)名。a和和b分別是定積分的下限和上限。分別是定積分的下限和上限。tol用來控制積分精度，缺省時取用來控制積分精度，缺省時取tol=0.001。trace控制是否展現(xiàn)積分過程，若取非控制是否展現(xiàn)積分過程，若取非0則展現(xiàn)積分過程，取則展現(xiàn)積分過程，取0則則不展現(xiàn)，缺省時取不展現(xiàn)，缺省時取trace=0。返回參數(shù)返回參數(shù)I即定積分值，即定積分值，n為被積函數(shù)的調(diào)用次數(shù)。為被積函數(shù)的調(diào)用次數(shù)。81MATLAB數(shù)值積分實例數(shù)值積分實例 (1) 建立被積函數(shù)文件建立被積函數(shù)文件fesin.m。function f=fesin(x)f=exp(-0.5*x).*sin(x+pi/6); (2) 調(diào)用數(shù)值積分函數(shù)調(diào)用數(shù)值積分函數(shù)quad求定積分。求定積分。S,n=quad(fesin,0,3*pi)S = 0.9008n = 7782quad8函數(shù)函數(shù)基于牛頓柯特斯法，