积分是非常重要的数学工具(微分方程、概率论等的基础),另一方

1、第九章 数值积分积分是非常重要的数学工具（微分方程、概率论等的基础），另一方面它们在实际问题中也有着许多直接的应用.函数的积分计算可分为数值积分和符号积分两类方法.数值积分是求积分的近似值的一类近似计算的方法.对于定积分d，如果对应的不定积分d不易求出，或者根本不能表示为初等函数时，那么就只能用数值积分的方法求其近似值了.例如，因为不定积分ed，d，ed，d无法由初等函数表示，所以计算这种类型的定积分只能用数值方法.至于由离散数据或图形表示的函数的定积分，理所当然地属于数值积分的范畴.符号积分是指求积分的精确值.本章重点介绍数值积分及其用MATLAB计算的方法.为了将数值积分的值与精确值比较，

2、在第一节中简单介绍用MATLAB 作符号积分的方法.9.1 积分的MATLAB符号计算计算机软件MATLAB的系统提供了符号计算定积分的函数int.下面主要介绍用MATLAB软件进行符号计算定积分、变上（下）限积分及其应用.9.1.1 定积分的MATLAB符号计算如果函数在区间上连续，且是的一个原函数，则函数在区间上的定积分为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 由，所围成的平面区域D.求平面区域D的面积S.解 （1） 画出由，所围成的封闭平面区域的图形，输入作函数图形的程序>> 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所围成的平面区域的图形')运行后屏幕显示图91.图91 由，所围成的平面区域图形（2）求平面区域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或数表，则常用的处理方法是利用各种插值函数（例如，样条插值函数、牛顿插值函数、拉格朗日插值多项式等）构造数值积分公式.本节主要介绍利用拉格朗日多项式和样条插值函数等构造数值积分公式的方法、梯形公式、辛普森（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)运行后屏幕显示取时，积分区间上等距节点的个数，用梯形公式计算定积分的值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系统提供了两种计算程序：一种是梯形数值积分（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主程序另外，也可以根据梯形公式和连续增加的

20、子区间数来逼近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 （复合辛普森公式）设函数在闭区间上的个等距节点处有定义，且其函数值为 ，则在上有, (9.11)使得 d, （9.12）其中是的截断误差，即余项.（9.11）式称为复合辛普森（Simpson）数值积分公式，简称辛普森（Simpson）公式.例9.3.7 用辛普森（Simpson）公式(9.1

25、1)计算ed，取个等距节点，并将计算结果与精确值比较，然后再取计算，观察对误差的影响.解 由，得.根据辛普森（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（取个等距节点）如下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.由此可见，节点的个数越大，步长越小，其绝对误差也越小.这个结论我们可以用证明推论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系统还提供了用自适应辛普森公式计算数值积分的程序quad.m，具体的调用格式

29、如下：调用格式一：quad(fun,a,b) 计算函数Y = fun(X)在区间(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，则会以动态的形式显示在递归求积分的整个过程的fcnt a b-a Q的值； 调用格式五：quad(fun,a,b, tol,TRACE,P1,P2, )此命令规定对函数fun (X,P1,P2,) 附加的参数P1, P2, .直接通过.传递tol 或 TRACE的空矩阵使用默认值.我们还可以根据复合辛普森（Simpson）公式编写名为comsimpson.m的计算d 的数值积分的MATLAB主程序.复合辛普森（Simpson）数值积分的MATLAB

31、主程序输入量：fun是被积函数， a,b分别是积分下、上限,是小区间的数目.输出量：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)在区间 (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，则会以动态的形式显示误差为tol的数值积分在递归求积分的整个过程和图形. 调用格式四：quad8(fun,a,b, tol,TRACE,P1,P2,.)此命令规定对函数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