数值积分及其MATLAB实现

9.1 积分的MATLAB符号计算9.1.1 定积分的MATLAB符号计算9.1.5 由，所围成的平面区域D.求平面区域D的面积S.解 输入作函数图形的程序 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), ylabel(y),title(y=sinx , y=cosx 和x=-0.5及x=1.5所围成的平面区域的图形)运行后屏幕显示图形.求平面区域D的面积S.输入计算面积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 及其近似值 Sj 如下S =2*2(1/2)+sin(1/2)-cos(1/2)-sin(3/2)-cos(3/2)Sj = 1.8269.1.2 变限积分的MATLAB符号计算例9.1.7 已知ed，求.解 输入程序： 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)运行后屏幕显示计算变限积分的导数.9.2 数值积分的思想及其MATLAB程序9.2.3 矩形公式的MATLAB程序（一） 函数sum的调用格式调用格式一：sum(X)调用格式二：sum (X,DIM)例9.2.2 用MATLAB和矩形公式（9.3）、（9.4）计算ed，并与精确值比较.解 将分成20等份，步长为，输入程序 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求定积分的精确值，输入程序 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）计算结果的绝对误差wz1、wz2.（二） 函数cumsum的调用格式调用格式一：cumsum(X) 调用格式二：cumsum (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(20)*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.3961wz1 = wz2 = 0.0088 0.00759.3 插值型数值积分及其MATLAB 程序9.3.2 梯形公式的MATLAB程序(一) 根据梯形公式和估计误差公式自己编写MATLAB程序计算定积分例9.3.2 分别取，用梯形公式计算定积分ed，并与精确值比较.然后观察对计算结果的有效数字和绝对误差的影响.解 编写并输入如下程序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(z8000-Fs)运行后屏幕显示取时，积分区间上等距节点的个数，用梯形公式计算定积分的值z8000和精确值intf的近似值Fs及其绝对误差Juewucha8000. (二) 用函数trapz计算定积分调用格式一：Z =trapz(Y)调用格式二：Z =trapz(X,Y)调用格式三：Z = trapz (X,Y,DIM) 或 trapz (Y,DIM) (三) 用函数cumtrapz计算定积分调用格式一：Z =cumtrapz (Y)调用格式二：Z =cumtrapz (X,Y)调用格式三：Z = cumtrapz (X,Y,DIM) 或 cumtrapz (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)*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.(四)梯形数值积分的MATLAB主程序梯形数值积分的MATLAB主程序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+h*s;endT=T(1:m);例9.3.6 用rctrap计算ed，递归14次，并将计算结果与精确值比较.解 输入程序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的绝对误差wT9.3.3 辛普森公式及其误差分析例9.3.7 用辛普森公式计算ed，取个等距节点，并将计算结果与精确值比较，然后再取计算，观察对误差的影响.解 由，得.根据辛普森（Simpson）公式编写并输入下面的程序 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(intf); Juewucha=abs(z-Fs)运行后屏幕显示用辛普森公式(9.11)计算定积分的近似值z和精确值intf及其绝对误差Juewucha（取个等距节点）.例9.3.8 估计用辛普森公式计算定积分ed时的误差，取.解 根据估计误差公式，先输入求的程序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*sin(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.2220e-0069.3.4 辛普森（Simpson）数值积分的MATLAB程序调用格式一：quad(fun,a,b)调用格式二：quad(fun,a,b,tol)调用格式三：Q,FCNT = quad (.) 调用格式四：quad(fun,a,b, tol,TRACE)调用格式五：quad(fun,a,b, tol,TRACE,P1,P2, )复合辛普森（Simpson）数值积分的MATLAB主程序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+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分别计算的近似值Q2、Q1和迭代次数FCNT14，取精度分别为，Q2、Q1分别与精确值Fs的绝对误差wQ2, wQ1如下9 0. 2.e-001 0. 11 0. 4.e-001 0. 13 0. 2.e-001 0.Q1 = FCNT14 = Q2 = 0.3413 13 0.3413Fs = wQ1 = wQ2 = 0.3413 3.6619e-009 3.7061e-0059.3.6 牛顿-科茨（Newton-Cotes）数值积分和误差分析的MATLAB程序(一) 估计误差的MATLAB程序计算n阶牛顿-科茨的公式的截断误差公式的MATLAB主程序function RNC=ncE(n)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);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= (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/80*(1/3*b-1/3*a)5*fxn1n = RNC4 = 4 8/945*(1/4*b-1/4*a)7*fxn2n = RNC8 = 8 543/88800*(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=diff(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)-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(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*myx8运行后屏幕显示误差限如下 RNC5 = RNC6 = 3.9320e-004 3.4823e-005(二) 计算科茨系数和求截断误差公式的MATLAB程序计算n阶科茨系数和求截断误差公式的MATLAB主程序function Cn, RNCn=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(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)*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);RNCn=(b-a)/n)(n+2)*fxn1*abs(y)/ suk;end例9.3.16 用计算n阶科茨系数和求截断误差公式的MATLAB主程序，计算定积分I=d的阶牛顿科茨公式的系数和截断误差公式.解 （1）先求阶牛顿科茨公式的系数和截断误差公式.输入程序 n1=1,Cn1,RNCn1=newcotE(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程序调用格式一：quad8(fun,a,b)调用格式二：quad8(fun,a,b,tol) 调用格式三： quad8(fun,a,b, tol,TRACE)调用格式四：quad8(fun,a,b, tol,TRACE,P1,P2,.)例9.3.17 用梯形公式、辛普森和牛顿科茨求积公式计算定积分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.742ztc = Columns 1 through 3 0 0.792 0.380.Columns 250 through 251 3.745 3.742（2）用辛普森求积公式计算定积分. 输入程序 syms xL= inline( exp(sin(x); QS,FCNTS =quad(L,0, pi/2,1.e-4,2)运行后屏幕显示用辛普森求积公式计算定积分的值QS和递归次数FCNTS分别如下QS = FCNTS = 3.254 13（3）用牛顿科茨求积公式计算定积分. 在MATLAB6.5中输入程序 syms xL= inline( exp(sin(x); Q8,FCNT8 = quad8(L,0, pi/2,1.e-4,3)运行后屏幕显示用牛顿科茨求积公式计算定积分的值Q8和递归次数FCNTS分别如下 Q8 = FCNT8 = 3.555 33（4）输入求定积分的精确值的程序 syms xy=exp(sin(x); F=int(y,0,pi/2), Fj=double(F)wzt=abs( Fj- zt), wQS = abs(Fj- QS), wQ8 = abs(Fj- Q8)运行后屏幕显示计算的定积分值F和近似值Fj，梯形公式、辛普森和牛顿科茨求积公式计算定积分的值与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.556 3.1430e-006wQS = wQ8 = 2.6399e-006 7.1127e-015（5）输入画图程序 syms xF=int(exp(sin(x),0,pi/2),Fj=double(F),plot(Fj,ro), hold onL= inline( exp(sin(x); Q=quad(L,0, pi/2,1.e-4,2),plot(Q,g*)hold off,hold on, h=pi/40;x=0:h:pi/2; y=exp(sin(x);ztc=cumtrapz(x,y), plot(x, ztc,mp), hold offhold on, Q8,FCNT8 = quad8(L,0, pi/2,1.e-4,3), hold offgrid, xlabel(自变量 X), ylabel(因变量 Y)title(牛顿科茨公式和梯形公式计算数值积分) legend(精确值, 辛普森公式计算定积分,梯形公式计算定积分,牛顿科茨公式计算定积分)运行后屏幕显示图形（略）.9.3.7 利用三次样条求表格型数值积分的MATLAB方法例9.3.18 给出概率积分e 的数据表： 0.46 0.47 0.48 0.49 f()0.322 9 0.319 9 0.316 8 0.313 8试首先用三次样条构造被积函数，再分别用quad和quadl函数计算定积分d，精度为，并将计算结果与精确值比较.解 方法1 首先用三次样条构造被积函数，再分别用quad和quadl函数计算定积分.在MATLAB工作窗口直接输入下面的程序 a=0.46,b=0.49, x =a: 0.01:b; y = 0.3229 0.3199 0.3168 0.3138, pp = spline(x,y) Q,FCNT= quad(ppval,a,b,pp)Ql,FCNTl= quadl(ppval,a,b,pp)运行后屏幕显示分别用辛普森和牛顿科茨公式计算的结果为pp = form: pp breaks: 0.4600 0.4700 0.4800 0.4900 coefs: 3x4 double pieces: 3 order: 4 dim: 1Q = FCNT = Ql = FCNTl = 0.0096 13 0.0096 18方法2 计算的精确解y及其近似解y1. 输入程序 syms xf=exp(-x.2)./( sqrt(2*pi); F= int(f, 0.46, 0.49)y=simple(F), y1=double(y)运行后屏幕显示结果y =42624/92261*pi(1/2)*(erf(49/100)-erf(23/50)y1 = 0.0096方法3 首先用三次样条构造被积函数，再分别用辛普森和牛顿科茨公式计算定积分.在MATLAB工作窗口输入程序X =0.46,0.47,0.48,0.49; Y= 0.3229 0.3199 0.3168 0.3138;s=interp1 (X,Y, X, spline); L3=poly2sym (s)运行后输出多项式L3. 保存名为L3.m的M文件function L3=L3(x)L3=3229/10000*x.3+3199/10000*x.2+198/625*x+1569/5000;然后输入程序Ql,FCNTl = quadl(L3, 0.46, 0.49,1.e-4) Q,FCNT = quad(L3, 0.46, 0.49,1.e-4)运行后屏幕显示计算结果Ql = FCNTl = Q = FCNT = 0.0171 18 0.0171 139.3.8 利用拉格朗日插值等方法求表格型数值积分的MATLAB方法例9.3.19 测得定积分d中被积函数在积分变量的某几个特定点X =0,0.157 1, 0.471 2, 0.549 8, 0.628 3, 0.863 9, 1.021 0, 1.099 6, 1.570 8处的值分别为Y= 0, 0.156 5, 0.454 0, 0.522 5, 0.587 8, 0.760 4, 0.852 6, 0.891 0, 1.000 0，试根据这些值首先分别用分段二、三阶拉格朗日插值多项式和三次样条插值函数构造被积函数，然后分别用quad函数和quad8函数计算，取精度为，并将计算结果与精确值1进行比较.解 方法1 (1) 拉格朗日插值多项式的MATLAB主程序function L=lfun(X,Y)m=length(X); n=M1; L=ones(m,m);for k=1:n+1 V=1; for i=1:n+1 if k=i V=conv(V,poly(X(i)/(X(k)-X(i); endendl(k,:)=poly2sym (V);endL=Y* l;(2) 用二阶拉格朗日插值多项式构造被积函数. 输入程序 X =0,0.1571, 0.4712,0.5498,0.6283,0.8639, 1.0210,1.0996,1.5708;Y=0,0.1565,0.4540,0.5225,0.5878,0.7604,0.8526,0.8910,1.0000;L1=lfun(X(1:3),Y(1:3), L2=lfun(X(3:5),Y(3:5),L3=lfun(X(5:7),Y(5:7), L4=lfun(X(7:9),Y(7:9),运行后输出多项式L1，L2，L3，L4.(3)输入程序 syms xL1=inline(-/*x.2+/*x);L2=inline(-/*x.2+/*x-/);L3=inline(-/*x.2+/*x-/);L4=inline(-/*x.2+/*x-/0); Q41,FCNT41 = quad(L1,0, 0.4712,1.e-4);Q42,FCNT42 =quad(L2, 0.4712, 0.6283,1.e-4);Q43,FCNT43 = quad(L3, 0.6283, 1.0210,1.e-4);Q44,FCNT44 = quad(L4, 1.0210, pi/2,1.e-4);Q= Q41+ Q42+ Q43+ Q44, FCNT= FCNT41+ FCNT42+ FCNT43+ FCNT44Q841,FCNT841 = quad8(L1,0, 0.4712,1.e-4);Q842,FCNT842 =quad8(L2, 0.4712, 0.6283,1.e-4);Q843,FCNT843 = quad8(L3, 0.6283, 1.0210,1.e-4);Q844,FCNT844 = quad8(L4, 1.0210, pi/2,1.e-4);Q8= Q841+ Q842+ Q843+ Q844, FCNT8= FCNT841+ FCNT842+ FCNT843+ FCNT844运行后屏幕显示分别用quad函数和quad8函数计算的近似值Q、Q8和迭代次数FCNT、FCNT8，取精度为如下Q = FCNT = Q8 = FCNT8 =0.413 52 0.413 132方法2 (1) 用三阶拉格朗日插值多项式构造被积函数.输入程序X=0,0.1571,0.4712,0.5498,0.6283,0.8639,1.0210,1.0996,1.5708;Y=0,0.1565,0.4540,0.5225,0.5878,0.7604,0.8526,0.8910,1.0000;L1=lfun(X(1:4),Y(1:4), L2=lfun(X(4:7),Y(4:7), L3=lfun(X(7:9),Y(7:9),运行后输出多项式L1，L2，L3. (2)输入程序 syms xL1=inline(-90323/*x.3-1389/*x.2+/*x);L2=inline(-/*x.3-42359/*x.2+/*x-/);L3 = inline(-/*x.2+/*x-/0); Q41,FCNT41 = quad(L1,0, 0.5498,1.e-4);Q42,FCNT42 =quad(L2, 0.5498, 1.0210,1.e-4);Q43,FCNT43 = quad(L3, 1.0210, pi/2,1.e-4);Q= Q41+ Q42+ Q43, FCNT= FCNT41+ FCNT42+ FCNT43Q841,FCNT841 = quad8(L1,0, 0.5498,1.e-4,2);Q842,FCNT842 =quad8(L2, 0.5498, 1.0210,1.e-4,2);Q843,FCNT843 = quad8(L3, 1.0210, pi/2,1.e-4,2);Q8= Q841+ Q842+ Q843, FCNT8= FCNT841+ FCNT842+ FCNT843运行后屏幕显示分别用quad函数和quad8函数计算的近似值Q、Q8和迭代次数FCNT、FCNT8，取精度为如下Q = FCNT= 0.584 39 Q84= FCNT8=0.584 99方法3 三次样条插值多项式计算数值积分.输入程序 a=0,b=1.5708, X=0,0.1571,0.4712,0.5498,0.6283,0.8639,1.0210,1.0996,1.5708;Y=0,0.1565,0.4540,0.5225,0.5878,0.7604,0.8526,0.8910,1.0000;pp = spline(X,Y) Q8,FCNT8= quad8(ppval,a,b, ,1.e-4,pp)Q,FCNT= quad(ppval,a,b,pp)运行后屏幕显示 Q8 = FCNT8 = Q = FCNT = 1.727 33 1.843 179.4 龙贝格（Romberg）公式及其MATLAB程序9.4.2 龙贝格积分的MATLAB程序龙贝格数值积分的MATLAB主程序function RT,R,wugu,h=romberg(fun,a,b, wucha,m)n=1;h=b-a; wugu=1; x=a;k=0; RT=zeros(4,4); RT(1,1)=h*(feval(fun,a)+feval(fun,b)/2;while(wuguwucha)&(km)|(kF=inline(1./(1+x); RT,R,wugu,h=romberg(F,0,1.5,1.e-8,13)syms x fi=int(1/(1+x),x,0,1.5); Fs=double(fi), wR=double(abs(fi-R), wR1= wR - wugu运行后屏幕显示精确值Fs，取精度为，利用龙贝格求积公式的romberg.m程序计算的近似值R，龙贝格积分表RT，误差估计err，最小步长h，近似值R与精确值Fs的绝对误差wR如下RT = Columns 1 through 4 1.000 0 0 00.143 0.857 0 00.828 0.057 0.470 00.571 0.152 0.292 0.400.717 0.766 0.540 0.2420.657 0.636 0.494 0.398 Columns 5 through 6 0 0 0 0 0 0 0 0 0.489 0 0.160 0.216R = wugu = h = 0.216 2.0648e-011 0.000Fs = wR = wR1 = 0.416 9.8125e-011 6.7478e-011如果用作为