積分是非常重要的數學工具(微分方程、概率論等的基礎),另一方

1、第九章 數值積分積分是非常重要的數學工具（微分方程、概率論等的基礎），另一方面它們在實際問題中也有著許多直接的應用.函數的積分計算可分為數值積分和符號積分兩類方法.數值積分是求積分的近似值的一類近似計算的方法.對于定積分d，如果對應的不定積分d不易求出，或者根本不能表示為初等函數時，那么就只能用數值積分的方法求其近似值了.例如，因為不定積分ed，d，ed，d無法由初等函數表示，所以計算這種類型的定積分只能用數值方法.至于由離散數據或圖形表示的函數的定積分，理所當然地屬于數值積分的范疇.符號積分是指求積分的精確值.本章重點介紹數值積分及其用MATLAB計算的方法.為了將數值積分的值與精確值比較，

2、在第一節(jié)中簡單介紹用MATLAB 作符號積分的方法.9.1 積分的MATLAB符號計算計算機軟件MATLAB的系統(tǒng)提供了符號計算定積分的函數int.下面主要介紹用MATLAB軟件進行符號計算定積分、變上（下）限積分及其應用.9.1.1 定積分的MATLAB符號計算如果函數在區(qū)間上連續(xù)，且是的一個原函數，則函數在區(qū)間上的定積分為d.我們可以用MATLAB軟件中的函數int對定積分進行符號計算，它的調用格式和功能如表91所示.表91 函數int求定積分的調用格式和功能函數int求定積分的調用格式功 能F=int(S,a,b)輸入量：S可以是被積函數，也可以是由幾個被積函數，構成的矩陣，積分上、下限

3、a 和 b可以是數，也可以是符號標量，積分變量是.輸出量：F是S關于積分變量的從a 到 b的定積分.F=int(S,t,a,b)輸入量：S、a 和 b同上，是積分變量.輸出量：F是S關于積分變量從a 到 b的定積分.例9.1.5 由，所圍成的平面區(qū)域D.求平面區(qū)域D的面積S.解 （1） 畫出由，所圍成的封閉平面區(qū)域的圖形，輸入作函數圖形的程序>> x=-1:0.001:2; F1= sin(x); F2=cos(x);plot(x ,F1,'b-',x ,F2,'g-'), axis(-1,pi/4+1,-1.3,1.3),xlabel('x

4、'), ylabel('y'),title('y=sinx , y=cosx 和x=-0.5及x=1.5所圍成的平面區(qū)域的圖形')運行后屏幕顯示圖91.圖91 由，所圍成的平面區(qū)域圖形（2）求平面區(qū)域D的面積S.解方程組得在內的兩條曲線交點的橫坐標，所求面積為dd.輸入計算面積S的程序 >> syms x f1= cos(x)-sin(x); f2=-f1; S1 =int(f1,x,-0.5,pi/4); S2=int(f2,x, pi/4,1.5); S=S1+S2,Sj= double (S)運行后屏幕顯示計算面積的值 S 及其近似值

5、Sj 如下S =2*2(1/2)+sin(1/2)-cos(1/2)-sin(3/2)-cos(3/2)Sj = 1.362037913188269.1.2 變限積分的MATLAB符號計算例9.1.7 已知ed，求.解 因為edteded所以輸入程序 >> syms x tF1=int(exp(t)*sin(2+sqrt(t3),x,0);F2=int(exp(t)*sin(2+sqrt(t3),0,x2);Fi= F1+ F2;dF=diff(Fi)運行后屏幕顯示計算變限積分的導數如下dF =-exp(x)*sin(2+(x3)(1/2)+2*x*exp(x2)*sin(2+(x

6、6)(1/2)9.2 數值積分的思想及其MATLAB程序9.2.3 矩形公式的MATLAB程序用矩形公式計算數值積分可用下面函數sum和cumsum實現，調用格式如下：（一） 函數sum的調用格式函數sum的兩種調用格式如下：調用格式一：sum(X)如果輸入向量X，則輸出向量X的所有元素的和，可用于按矩形公式（9.3）和（9.4）計算積分；如果輸入矩陣X，則輸出矩陣X的每列向量的和.調用格式二：sum (X,DIM)輸入N-D 數組X，輸出為X的每個DIM向量的和.例9.2.2 用MATLAB和矩形公式（9.3）、（9.4）計算ed，并與精確值比較.解 將分成20等分，步長為，輸入程序>

7、> h=pi/40; x=0:h:pi/2; y=exp(sin(x);z1=sum(y(1:20)*h, z2=sum(y(2:21)*h,運行后屏幕顯示矩形公式計算結果分別如下 z1 = z2 = 3.0364 3.1713即矩形公式（9.3）計算結果3.036 4, 矩形公式（9.4）計算結3.171 3.求定積分的精確值，輸入程序>> syms x F=int(exp(sin(x),x,0, pi/2), Fs= double (F), wz1=abs( Fs-z1), wz2= abs( Fs-z2)運行后屏幕顯示定積分的精確值Fs和與用矩形公式（9.3）、（9.4

8、）計算結果的絕對誤差wz1、wz2分別如下Warning: Explicit integral could not be found.> In C:MATLAB6p5p1toolboxsymbolicsymint.m at line 58F = Fs =int(exp(sin(x),x = 0 . 1/2*pi) 3.1044wz1 = wz2 = 0.0680 0.0670由此可見，該定積分用矩形公式（9.4）計算結果比矩形公式（9.3）的絕對誤差小.（二） 函數cumsum的調用格式函數cumsum的調用格式如下：調用格式一：cumsum(X) 輸入向量X，輸出向量X的元素依次累加和

9、，可用于按矩形公式（9.3）、（9.4）計算積分； 輸入矩陣X，輸出矩陣為X的每列向量的依次累加和.調用格式二：cumsum (X,DIM)輸入N-D 數組X，輸出為X的每個DIM向量的依次累加和.例9.2.4 用MATLAB的函數sum 和 cumsum及矩形公式（9.3）、（9.4）計算ed，并與精確值比較.解 將分成20等分，步長為，輸入程序如下（注意sum 和 cumsum的用法）>> h=pi/40; x=0:h:pi/2; y=exp(-x).*sin(x); z1=sum(y(1:20)*h,z2=sum(y(2:21)*h, z=cumsum(y); z11=z(2

10、0)*h, z12=(z(21)-z(1)*h,運行后屏幕顯示計算結果分別如下z1 = z2 = z11 = z12 = 0.3873 0.4036 0.3873 0.4036求定積分的精確值，輸入程序>> syms x F=int(exp(-x)*sin(x),x,0, pi/2)Fs= double (F) ,wz1=abs( Fs-z1), wz2= abs( Fs-z2)運行后屏幕顯示定積分的精確值Fs和用矩形公式（9.3）、（9.4）計算結果的絕對誤差wz1、wz2分別如下F = Fs =1/2*(-1+exp(pi)(1/2)/exp(pi)(1/2) 0.3961wz

11、1 = wz2 = 0.0088 0.0075由此可見，用矩形公式（9.4）計算該定積分的結果比用矩形公式（9.3）的結果的絕對誤差小.9.3 插值型數值積分及其MATLAB 程序如果被積函數和積分變量是離散的數據X，Y或數表，則常用的處理方法是利用各種插值函數（例如，樣條插值函數、牛頓插值函數、拉格朗日插值多項式等）構造數值積分公式.本節(jié)主要介紹利用拉格朗日多項式和樣條插值函數等構造數值積分公式的方法、梯形公式、辛普森（Simpson）公式和一般的Newton-Cotes公式及其誤差分析，三者之間的關系，并通過實例說明用MATLAB軟件計算數值積分的方法.9.3.2 梯形公式的MATLAB程

12、序(一) 根據梯形公式和估計誤差公式自己編寫MATLAB程序計算定積分例9.3.2 分別取，用梯形公式計算定積分ed，并與精確值比較.然后觀察對計算結果的有效數字和絕對誤差的影響.解 根據（9.5a）式編寫并輸入如下程序>> h=pi/8000; a=0;b=pi/2; x=a:h:b;n=length(x), y=exp(sin(x);z1=(y(1)+y(n)*h/2; z2=sum(y(2:n-1)*h; z8000=z1+z2,syms tf=exp(sin(t); intf=int(f,t,a,b), Fs=double(intf),Juewucha8000=abs(z8

13、000-Fs)運行后屏幕顯示取時，積分區(qū)間上等距節(jié)點的個數，用梯形公式計算定積分的值z8000和精確值intf的近似值Fs及其絕對誤差Juewucha8000如下n = z8000 = 4001 3.10437900500451Warning: Explicit integral could not be found.> In C:MATLAB6p5toolboxsymbolicsymint.m at line 58intf = Fs =int(exp(sin(t),t = 0 . 1/2*pi) 3.10437901785556Juewucha8000 = 1.285104200832

14、166e-008然后再分別取，運行后屏幕依次顯示用梯形公式計算定積分的值z800、z80和與精確值intf的近似值Fs的絕對誤差Juewucha 800、Juewucha 80如下n800 = z800 = Juewucha800 = 401 3.10437773275082 1.285104739068288e-006n80 = z80= Juewucha80 =41 3.10425050738255 1.285104730022191e-004 (二) 用函數trapz計算定積分用MATLAB軟件和梯形公式計算數值積分和估計誤差限除了用例9.3.2和例9.3.3自己編寫的MATLAB計算程

15、序外，MATLAB系統(tǒng)提供了兩種計算程序：一種是梯形數值積分（Trapezoidal doubleal integration.）的程序trapz.m；另一種是累加梯形數值積分（Cumulative trapezoidal doubleal integration.）的程序cumtrapz.m，它們的調用格式如下：調用格式一：Z =trapz(Y)輸入Y，輸出為按梯形公式（9.5a）計算的Y的積分的近似值（單位步長）.（1）如果輸入Y是向量，則Z是Y的積分；（2）如果輸入Y是矩陣，則行向量Z是Y的每列的積分.調用格式二：Z =trapz(X,Y)輸入X,Y為同長度的數組，輸出為按梯形公式（9.

16、5a）計算Y對X的積分（步長不一定相等）.調用格式三：Z = trapz (X,Y,DIM) 或 trapz (Y,DIM)輸出為按梯形公式（9.5a）計算Y 的每個DIM的積分，其中X的長度必須和size(Y,DIM)相同.(三) 用函數cumtrapz計算定積分調用格式一：Z =cumtrapz (Y)輸入Y，輸出用梯形方法計算的Y的累加積分（Cumulative integral）的近似值（單位步長）.（1）如果輸入Y是向量，則Z是Y的累加積分；（2）如果輸入Y是矩陣，則行向量Z是Y的每列的累加積分.調用格式二：Z =cumtrapz (X,Y)輸入X,Y為同長度的數組，輸出用梯形方法計

17、算的Y對X的累加積分（步長不一定相等）.調用格式三：Z = cumtrapz (X,Y,DIM) 或 cumtrapz (Y,DIM)輸出用梯形方法計算的Y 的每個DIM的累加積分，其中X的長度必須和size(Y,DIM)相同. 例9.3.4 用MATLAB的函數trapz 和cumtrapz分別計算ed，精確到，并與矩形公式（9.3），（9.4）比較.解 將分成20等分，步長為，輸入程序如下（注意trapz(y)是單位步長, trapz(y)*h=trapz(x,y)）：>> h=pi/40; x=0:h:pi/2; y=exp(-x).*sin(x);z1=sum(y(1:20

18、)*h, z2=sum(y(2:21)*h, z=(z1+z2)/2z3=trapz(y)*h, z3h=trapz(x,y), z3c=cumtrapz(y)*h,運行后屏幕顯示用矩形公式（9.3），（9.4）計算結果z1、z2和二者的平均數z、函數trapz 和cumtrapz分別計算結果z3、z3c分別如下z1 = z2 = z = z3 = z3h =0.3873 0.4036 0.3954 0.3954 0.3954z3c = Columns 1 through 7 0 0.0028 0.0109 0.0234 0.0395 0.0586 0.0798 Columns 8 throu

19、gh 14 0.1028 0.1270 0.1519 0.1771 0.2023 0.2273 0.2517 Columns 15 through 21 0.2755 0.2983 0.3201 0.3408 0.3602 0.3785 0.3954顯然，函數trapz 和cumtrapz分別計算結果都與矩形公式（9.3），（9.4）計算結果的平均數z相等，即z3=z3c=z=0.395 4，梯形公式計算的結果與精確值Fs = 0.396 1更接近.但是，函數cumtrapz計算的結果展示了每次累加計算過程的結果.(四) 自編梯形數值積分的MATLAB主程序另外，也可以根據梯形公式和連續(xù)增加的

20、子區(qū)間數來逼近d 第次遞歸在個等距點處對的值的和的方法，現提供計算數值積分的程序如下：梯形數值積分的MATLAB主程序輸入量：fun是被積函數， a,b分別是積分下、上限，是遞歸的次數.輸出量：T是利用梯形公式數值計算d的遞歸值.function T=rctrap(fun,a,b,m)n=1;h=b-a; T=zeros(1,m+1); x=a; T(1)=h*(feval(fun,a)+feval(fun,b)/2;for i=1:m h=h/2; n=2*n; s=0; for k=1:n/2 x=a+h*(2*k-1); s=s+feval(fun,x);endT(i+1)=T(i)/2

21、+h*s;endT=T(1:m);例9.3.6 用rctrap計算ed，遞歸14次，并將計算結果與精確值比較.解 （1）建立fun.m文件命名的M文件函數function y = fun(x) y = exp( (-x.2)./2)./(sqrt(2*pi);（2）輸入程序>>T=rctrap(fun,0,pi/2,14), syms t fi=int(exp(-t2)/2)/(sqrt(2*pi),t,0, pi/2); Fs= double(fi), wT= double(abs(fi-T)運行后屏幕顯示精確值Fs，用rctrap計算的遞歸值T和T與精確值Fs的絕對誤差wT如下

22、T = Columns 1 through 4 0.40457385587044 0.43245899179539 0.43953669400491 0.44129851884975 Columns 5 through 8 0.44173842995173 0.44184837274713 0.44187585624611 0.44188272698315 Columns 9 through 12 0.44188444465881 0.44188487407718 0.44188498143174 0.44188500827038 Columns 13 through 14 0.4418850

23、1498004 0.44188501665745Fs = 0.44188501721659wT = Columns 1 through 4 0.03731116134615 0.00942602542121 0.00234832321168 0.00058649836684 Columns 5 through 8 0.00014658726486 0.00003664446946 0.00000916097048 0.00000229023344 Columns 9 through 12 0.00000057255779 0.00000014313941 0.00000003578485 0.

24、00000000894621 Columns 13 through 14 0.00000000223655 0.00000000055914由此可見，遞歸14次，用rctrap計算數值積分的結果與精確值的絕對誤差為，0.441 885 017.9.3.3 辛普森（Simpson）公式及其誤差分析定理9.2 （復合辛普森公式）設函數在閉區(qū)間上的個等距節(jié)點處有定義，且其函數值為 ，則在上有, (9.11)使得 d, （9.12）其中是的截斷誤差，即余項.（9.11）式稱為復合辛普森（Simpson）數值積分公式，簡稱辛普森（Simpson）公式.例9.3.7 用辛普森（Simpson）公式(9.1

25、1)計算ed，取個等距節(jié)點，并將計算結果與精確值比較，然后再取計算，觀察對誤差的影響.解 由，得.根據辛普森（Simpson）公式(9.11)編寫并輸入下面的程序>> a=0;b=1;m=10000; h=(b-a)/(2*m); x=a:h:b; y=exp(-x.2)./2)./(sqrt(2*pi);z1=y(1)+y(2*m+1); z2=2*sum(y(2:2:2*m); z3=4*sum(y(3:2:2*m);z=(z1+z2+z3)*h/3, syms t,f=exp(-t2)/2)/(sqrt(2*pi);intf=int(f,t,a,b), Fs=double(i

26、ntf); Juewucha=abs(z-Fs)運行后屏幕顯示用辛普森公式(9.11)計算定積分的近似值z和精確值intf及其絕對誤差Juewucha（取個等距節(jié)點）如下z = 0.34133406408431intf =1125899906842624/5644425081792261*erf(1/2*2(1/2)*2(1/2)*pi(1/2)Juewucha = 1.068198423692657e-005然后再取計算的近似值z和與精確值intf及其絕對誤差Juewucha的結果為：z = 0.323 261 353 494 30, Fs = 0.341 344 746 068 54, J

27、uewucha = 0.018 083 392 574 25.由此可見，節(jié)點的個數越大，步長越小，其絕對誤差也越小.這個結論我們可以用證明推論9.1的方法證明，其結果見下面的推論.例9.3.8 估計用辛普森公式計算定積分ed時的誤差，取.解 根據估計誤差公式（9.14），先輸入求的程序>>syms x,y=exp(sin(x); yx4=diff(y,x,4)運行后輸出被積函數的四階導函數. 然后在輸入誤差估計程序>>h=pi/40; x=0:0.00001:pi/2; yx4=sin(x).*exp(sin(x)-4*cos(x).2.*exp(sin(x)+3*si

28、n(x).2.*exp(sin(x)-6*sin(x).*cos(x).2.*exp(sin(x)+cos(x).4.*exp(sin(x);juyx4= abs(yx4); RS=(h4)*(pi/2)*max(juyx4)/180運行后屏幕顯示誤差估計值RS = 3.610450295892220e-0069.3.4 辛普森（Simpson）數值積分的MATLAB程序用MATLAB軟件和辛普森（Simpson）公式計算數值積分和估計誤差限除了用例9.3.7和例9.3.8自己編寫的MATLAB計算程序外，MATLAB系統(tǒng)還提供了用自適應辛普森公式計算數值積分的程序quad.m，具體的調用格式

29、如下：調用格式一：quad(fun,a,b) 計算函數Y = fun(X)在區(qū)間(a,b)上的數值積分，自動選擇步長，在MATLAB 6.3中絕對誤差限為，在MATLAB 5.3中誤差限為，輸出數值積分值.其中函數Y = fun(X)是以fun.m文件命名的M文件函數（或庫函數如sin,log），或者是F = inline(fun)形式，或者數組，函數Y = fun(X)將接受向量的自變量，并且返回結果是向量Y，即在X的每個元素處的積分估計值.調用格式二：quad(fun,a,b,tol) 同上，但指定絕對誤差限為tol，默認值為 ，在MATLAB 5.3中誤差限為.調用格式三：Q,FCNT

30、= quad (.) 輸出數值積分值Q和函數賦值的編號FCNT.調用格式四：quad(fun,a,b, tol,TRACE)輸入量非零數TRACE，則會以動態(tài)的形式顯示在遞歸求積分的整個過程的fcnt a b-a Q的值； 調用格式五：quad(fun,a,b, tol,TRACE,P1,P2, )此命令規(guī)定對函數fun (X,P1,P2,) 附加的參數P1, P2, .直接通過.傳遞tol 或 TRACE的空矩陣使用默認值.我們還可以根據復合辛普森（Simpson）公式編寫名為comsimpson.m的計算d 的數值積分的MATLAB主程序.復合辛普森（Simpson）數值積分的MATLAB

31、主程序輸入量：fun是被積函數， a,b分別是積分下、上限,是小區(qū)間的數目.輸出量：y是利用復合辛普森數值積分公式(9.11)計算d的數值積分結果.function y=comsimpson(fun,a,b,n) z1=feval (fun,a)+ feval (fun,b);m=n/2;h=(b-a)/(2*m); x=a;z2=0; z3=0; x2=0; x3=0;for k=2:2:2*m x2=x+k*h; z2= z2+2*feval (fun,x2); endfor k=3:2:2*m x3=x+k*h; z3= z3+4*feval (fun,x3); endy=(z1+z2+

32、z3)*h/3;例9.3.9 用comsimpson.m和quad.m分別計算定積分ed，取精度為，并與精確值比較.解 輸入程序>> Q1,FCNT14 = quad(fun,0,1,1.e-4,3), Q2 =comsimpson (fun,0,1,10000)syms x fi=int(exp( (-x.2)./2)./(sqrt(2*pi),x,0, 1); Fs= double (fi)wQ1= double (abs(fi-Q1) ), wQ2= double (abs(fi-Q2) )運行后屏幕顯示的精確值Fs，用comsimpson.m和quad.m分別計算的近似值Q

33、2、Q1和迭代次數FCNT14，取精度分別為，Q2、Q1分別與精確值Fs的絕對誤差wQ2, wQ1如下9 0.0000000000 2.71580000e-001 0.1070275100 11 0.2715800000 4.56840000e-001 0.1597942242 13 0.7284200000 2.71580000e-001 0.0745230082Q1 = FCNT14 = Q2 = 0.3413 13 0.3413Fs = wQ1 = wQ2 = 0.3413 3.6619e-009 3.7061e-0059.3.6 牛頓-科茨（Newton-Cotes）數值積分和誤差分析

34、的MATLAB程序下面介紹計算科茨（Cotes）系數和估計（9.21）式在上的截斷誤差的程序及其用法，并且介紹根據在MATLAB中提供了用自適應遞歸牛頓科茨公式計算數值積分的程序quad8.m和使用方法.(一) 估計誤差的MATLAB程序根據計算定積分d的近似值的n階牛頓-科茨公式的截斷誤差公式（9.24）式和（9.25）式,現提供的求（9.21）式在上的截斷誤差公式的主程序如下：計算n階牛頓-科茨的公式的截斷誤差公式的MATLAB主程序輸入量： 牛頓-科茨的公式（9.21）的階數.輸出量： RNC是對應的n階牛頓-科茨的公式（9.21）的截斷誤差公式.function RNC=ncE(n)s

35、uk=1; p=n/2-fix(n/2);if p=0for k=1:n+2suk=suk*k;endsuk; syms t a b fxn2,su=t2; for u=1:nsu=su*(t-u);endsu; intf=int(su,t,0,n); y=double(intf);RNC= (b-a)/n)(n+3)*fxn2*abs(y)/ suk;elsefor k=1:n+1suk=suk*k;endsuk; syms t a b fxn1,su=t; for u=1:nsu=su*(t-u);endsu; intf=int(su,t,0,n); y=double(intf);RNC=

36、 (b-a)/n)(n+2)*fxn1*abs(y)/ suk;end例9.3.14 用求截斷誤差公式的MATLAB主程序，求計算定積分d的近似值的階牛頓科茨公式的截斷誤差公式.解 輸入求階牛頓科茨公式的截斷誤差公式的程序>> n=1, RNC1=ncE(n), n=2, RNC2=ncE(n), n=3, RNC3=ncE(n)n=4, RNC4=ncE(n), n=8, RNC8=ncE(n)運行后屏幕顯示結果如下 n = RNC1 = 1 1/12*(b-a)3*fxn1n = RNC2 = 2 1/90*(1/2*b-1/2*a)5*fxn2n = RNC3 = 3 3/8

37、0*(1/3*b-1/3*a)5*fxn1n = RNC4 = 4 8/945*(1/4*b-1/4*a)7*fxn2n = RNC8 = 8 35065906189543/6926923254988800*(1/8*b-1/8*a)11*fxn2例9.3.15 用MATLAB軟件估計用5、6階牛頓科茨公式計算定積分ed的截斷誤差公式和誤差限，取15位小數.解 輸入求階牛頓科茨公式的截斷誤差公式和被積函數的6，8階導函數的程序>> n=5;RNC5=ncE(n), n=6;RNC6=ncE(n)syms xy=exp(sin(x); yx6=diff(y,x,6), yx8=dif

38、f(y,x,8)運行后輸出被積函數的6，8階導函數（略）和階牛頓科茨公式的截斷誤差公式如下 RNC5 = RNC6 =275/12096*(1/5*b-1/5*a)7*fxn1 9/1400*(1/6*b-1/6*a)9*fxn2然后再輸入誤差估計程序>>a=0;b=pi/2; h=pi/40; x=0:0.00001:pi/2;yx6=-sin(x).*exp(sin(x)+16*cos(x).2.*exp(sin(x)-15*sin(x).2.*exp(sin(x)+75*sin(x).*cos(x).2.*exp(sin(x)-20*cos(x).4.*exp(sin(x)-

39、15*sin(x).3.*exp(sin(x)+45*sin(x).2.*cos(x).2.*exp(sin(x)-15*sin(x).*cos(x).4.*exp(sin(x)+cos(x).6.*exp(sin(x);yx8=cos(x).8.*exp(sin(x)-756*sin(x).*cos(x).2.*exp(sin(x)-1260*sin(x).2.*cos(x).2.*exp(sin(x)+630*sin(x).*cos(x).4.*exp(sin(x)-420*sin(x).3.*cos(x).2.*exp(sin(x)+210*sin(x).2.*cos(x).4.*exp

40、(sin(x)-28*sin(x).*cos(x).6.*exp(sin(x)-56*cos(x).6.*exp(sin(x)+sin(x).*exp(sin(x)+63*sin(x).2.*exp(sin(x)+210*sin(x).3.*exp(sin(x)+105*sin(x).4.*exp(sin(x)-64*cos(x).2.*exp(sin(x)+336*cos(x).4.*exp(sin(x);myx6=max(yx6); myx8=max(yx8); RNC5 =275/12096*(1/5*b-1/5*a)7*myx6RNC6 =9/1400*(1/6*b-1/6*a)9*m

41、yx8運行后屏幕顯示誤差限如下 RNC5 = RNC6 = 3.625385713339320e-004 3.826183594914823e-005(二) 計算科茨（Cotes）系數和求截斷誤差公式的MATLAB程序現提供計算定積分I=d的n階牛頓科茨的公式（9.21）、（9.23）編寫的計算科茨（Cotes）系數主程序如下：計算n階科茨（Cotes）系數和求截斷誤差公式的MATLAB主程序輸入量：牛頓科茨的公式（9.21）的階數.輸出量：Cn是計算定積分I的n階牛頓科茨的公式的科茨（Cotes）系數，RNCn是對應的n階牛頓科茨的公式（9.21）的截斷誤差公式.function Cn, R

42、NCn=newcotE(n)syms t a b M,Fz=zeros(1,n+1); Cn=zeros(1,n+1); su=t;k=1;m=1;m0=1; for u=1:nsu1=su*(t-u);m01=m0*u; su=su1;m0=m01;endsu;m0; f1=su/(t-0); intf1=int(f1,t,0,n); y= double(intf1);Cn(1)=(-1)(n-0)*y)/(n*m0); k=1;m=1;for j=1:nk1=k*j; m1=m*(n-j); f=su/(t-j); intf=int(f,t,0,n); y=double(intf);Cn(

43、j+1)=(-1)(n-j)*y)/(n*k1*m1);warning off MATLAB:divideByZeroend fn=su/(t-n); intfn=int(fn,t,0,n); y= double(intfn);Cn(n+1)=y/(n*m0);Cn; suk=1; p=n/2-fix(n/2);if p=0for k=1:n+2suk=suk*k;endsuk; syms t a b fxn2,su=t2; for u=1:nsu=su*(t-u);endsu; intf=int(su,t,0,n); y=double(intf);RNCn=(b-a)/n)(n+3)*fxn

44、2*abs(y)/suk;elsefor k=1:n+1suk=suk*k;endsuk; syms t a b fxn1,su=t; for u=1:nsu=su*(t-u);endsu; intf=int(su,t,0,n); y=double(intf);RNCn=(b-a)/n)(n+2)*fxn1*abs(y)/ suk;end例9.3.16 用計算n階科茨（Cotes）系數和求截斷誤差公式的MATLAB主程序，計算定積分I=d的階牛頓科茨公式的系數和截斷誤差公式.解 （1）先求階牛頓科茨公式的系數和截斷誤差公式.輸入程序>> n1=1,Cn1,RNCn1=newcotE

45、(n1)n2=2,Cn2,RNCn2=newcotE(n2)n3=3,Cn3,RNCn3=newcotE(n3)運行后屏幕顯示階牛頓科茨公式的系數Cn1, Cn2, Cn3和截斷誤差公式RNCn1, RNCn2, RNCn3如下n1 = 1 Cn1 = 1/2 1/2 RNCn1 =1/12*(b-a)3*fxn1n2 = 2 Cn2 = 1/6 2/3 1/6 RNCn2 =1/90*(1/2*b-1/2*a)5*fxn2n3 = 3 Cn3 = 1/8 3/8 3/8 1/8 RNCn3 =3/80*(1/3*b-1/3*a)5*fxn1(三) 計算牛頓科茨公式的MATLAB程序在MATL

46、AB5.3中提供了用自適應遞歸牛頓科茨公式計算數值積分的程序quad8.m.但是，因為quad8.m遞歸速度慢和其它原因，在MATLAB6.5和MATLAB7.01中極力推薦使用程序quadl.m取代quad8.m. 在MATLAB7.01中雖然還保留quad8.m的文件名，實際上已經將它的程序廢棄，輸入有關quad8的程序時，運行的是quadl.m的程序.quad8.m和quadl.m的調用方法與quad 類似，下面介紹quad8.m的調用格式.調用格式一：quad8(fun,a,b) 計算函數Y = fun(X)在區(qū)間 (a,b)上的數值積分，自動選擇步長，相對誤差限為，輸出數值積分值.函

47、數Y = fun(X)將接受向量的自變量X，并且返回結果是向量Y，即在X的每個元素處的積分估計值.調用格式二：quad8(fun,a,b,tol) 同上，返回相對誤差tol的數值積分，如果TOL = rel_tol abs_tol，則rel_tol 和abs_tol 分別表示相對誤差和絕對誤差.調用格式三： quad8(fun,a,b, tol,TRACE)輸入非零數TRACE，則會以動態(tài)的形式顯示誤差為tol的數值積分在遞歸求積分的整個過程和圖形. 調用格式四：quad8(fun,a,b, tol,TRACE,P1,P2,.)此命令規(guī)定對函數fun (X,P1,P2,.) 附加的參數P1,

48、P2, .直接通過.傳遞tol 或 TRACE的空矩陣使用默認值.例9.3.17 用梯形公式、辛普森（Simpson）和牛頓科茨求積公式計算定積分ed，取精度，作出它們的積分圖，并與精確值進行比較;解 （1）用梯形求積公式計算定積分. 輸入程序>> h=pi/500; x=0:h:pi/2; y=exp(sin(x);zt=trapz(x,y), ztc=cumtrapz(x,y), plot(x, ztc,'ro')運行后屏幕顯示用函數trapz 和cumtrapz分別計算結果zt、ztc分別如下zt = 3.10437572798742ztc = Columns

49、 1 through 3 0 0.00630298652792 0.01264569951380.Columns 250 through 251 3.08729642810745 3.10437572798742（2）用辛普森（Simpson）求積公式計算定積分. 輸入程序>> syms xL= inline(' exp(sin(x)'); QS,FCNTS =quad(L,0, pi/2,1.e-4,2)運行后屏幕顯示用辛普森（Simpson）求積公式計算定積分的值QS和遞歸次數FCNTS分別如下QS = FCNTS = 3.10438133817254 13（3

50、）用牛頓科茨求積公式計算定積分. 在MATLAB6.5中輸入程序>> syms xL= inline(' exp(sin(x)'); Q8,FCNT8 = quad8(L,0, pi/2,1.e-4,3)運行后屏幕顯示用牛頓科茨求積公式計算定積分的值Q8和遞歸次數FCNTS分別如下 Q8 = FCNT8 = 3.10437901785555 33（4）輸入求定積分的精確值的程序>> syms xy=exp(sin(x); F=int(y,0,pi/2), Fj=double(F)wzt=abs( Fj- zt), wQS = abs(Fj- QS), w

51、Q8 = abs(Fj- Q8)運行后屏幕顯示計算的定積分值F和近似值Fj，梯形公式、辛普森（Simpson）和牛頓科茨求積公式計算定積分的值與Fj的絕對誤差wztc, wQS和wQ8如下Warning: Explicit integral could not be found.> In C:MATLAB6p5p1toolboxsymbolicsymint.m at line 58F =int(exp(sin(x),x = 0 . 1/2*pi)Fj = wzt = 3.10437901785556 3.289868133471430e-006wQS = wQ8 = 2.32031698