第四章数值积分第九章

插值型数值积分及其程梯 (一)根据梯 和估计误 程序计算定积

/例4.1.1分别取h/8000,/800,/80，用梯 计算定积分I 解>>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,symsf=exp(sin(t));intf=int(f,t,a,b),Fs=double(intf),运行后屏幕显示取h8000时，积分区间[02上等距节点的个数nIz8000intfFstrapz计算定积分调用格式一：Z=trapz(Y)调用格式二：Z调用格式三：ZtrapzX,Y,DIM)trapzcumtrapz计算定积分调用格式一：Z=cumtrapz(Y)调用格式二：ZcumtrapzX,Y)调用格式三：ZcumtrapzX,Y,DIM)cumtrapz/例4.1.2 的函数trapz和cumtrapz分别计算

exsinxdx解将[02]20等份，步长为40，输入程序如下（trapz(y)是单位步长,trapz(y)*h=trapz(x,y)：>>h=pi/40;x=0:h:pi/2;y=exp(-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,运行后屏幕显示用矩形计算结果z1z2和二者的平均数z函数trapz和 functionT=rctrap(fun,a,b,m)n=1;h=b-a;T=zeros(1,m+1);x=a;fori=1:mh=h/2;n=2*n;s=0;fork=1:n/2x=a+h*(2*k-1);s=s+feval(fun,x);例4.1.3用rctrap计算I 值比较解

2e2dx140>>T=rctrap(@fun,0,pi/2,14),symstfi=int(exp((-t^2)/2)/(sqrt(2*pi)),t,0,pi/2);Fs=double(fi),wT=double(abs(fi-T))111例4.1.4 计算I 0

e2dx，取n20001计算结果与精确值比较，然后再取n13计算，观察n对误差的影响）解由n2m120001，得m10000.根据 （Simpson编写并输入）>>a=0;b=1;m=10000;h=(b-a)/(2*m);z1=y(1)+y(2*m+1);z2=2*sum(y(2:2:2*m));z=(z1+z2+z3)*h/3,symst,f=exp((-t^2)/2)/(sqrt(2*pi));intf=int(f,t,a,b),Fs=double(intf);Juewucha=abs(z-运行后屏幕显示用(9.11)计算定积分I的近似值z和精确值intf及其绝对误Juewucha（取n20001个等距节点例4.1.5估计用计算定积分I2

esinxdx时的误差，取h40解根据估计误差，先输入求f(4)(x)的程>>symsx,y=exp(sin(x));运行后输出被积函数的四阶导函数.>>h=pi/40;juyx4=abs(yx4);RS=(h^4)*(pi/2)*max(juyx4)/180RS调用格式一调用格式二调用格式三 T]=quad调用格式五：quad(‘fun’,a,b,tol,TRACE,P1,P2,复合（Simpson）数值积分的主程 z1=feval(fun,a)+feval(fun,b);m=n/2;h=(b-a)/(2*m);x=a;z2=0;z3=0;x2=0;forx2=x+k*h;z2=z2+2*fevalforx3=x+k*h;z3=z3+4*feval(fun,x3);例4.1.6用comsimpson.m和quad.m分别计算定积分I 解>> T14]=quad(@fun,0,1,1.e-Q2=comsimpson(@fun,0,1,10000)symsx

1e10

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)IFscomsimpson.mquad.mIQ2、Q1FCNT14，取精度分别为104，Q2、Q1FswQ2,wQ1如下9e-e-Q1FCNT14Q2=FswQ1wQ2 (一)估计误差 程 的截断误差的主程functionRNC=ncE(n)suk=1;p=n/2-fix(n/2);ifp==0fork=1:n+2suk;symstabfxn2,su=t^2;foru=1:nsu;intf=int(su,t,0,n);y=double(intf);RNC=(((b-a)/n)^(n+3))*fxn2*abs(y)/suk;fork=1:n+1suk;symstabfxn1,su=t;foru=1:nsu;intf=int(su,t,0,n);y=double(intf);RNC=(((b-a)/n)^(n+2))*fxn1*abs(y)/suk;b例4.1.7用求截断误差 主程序，求计算定积分af(x)dx的近似值的1,2,3,4,8阶Newton-Cotes的截断误差.b>>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)nn1RNC1RNC2=n23RNC3=nn48RNC4RNC8=例4.1.8用软件估计用5、6阶Newton-Cotes计算定积分I20解输入求n5,6阶Newton-Cotes的截断误差和被积函数的6，8阶导>>n=5;RNC5=ncE(n),symsy=exp(sin(x));yx6=diff(y,x,6),运行后输出被积函数的6，8阶导函数（略）和n5,6阶Newton-Cotes的截断误RNC5 RNC6 >>a=0;b=pi/2;h=pi/40;myx6=max(yx6);myx8=max(yx8);RNC5=275/12096*(1/5*b-RNC6=9/1400*(1/6*b-RNC5 RNC6 k(二)计算cotes系数C(n)和求截断误差的程kk计算n阶cotes系数C(n)和求截断误差的主程kfunction[Cn, symstabM,Fz=zeros(1,n+1);Cn=zeros(1,n+1);su=t;k=1;m=1;m0=1;forsu1=su*(t-u);m01=m0*u;su=su1;m0=m01;su;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;forj=1:nk1=k*j;m1=m*(n-j);f=su/(t-j);intf=int(f,t,0,n);y=double(intf);warningoff:divideByZerofn=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);iffork=1:n+2suk;symstabfxn2,su=t^2;foru=1:nsu;intf=int(su,t,0,n);fork=1:n+1suk;symstabfxn1,su=t;foru=1:nsu;intf=int(su,t,0,n);=(((b-a)/n)^(n+2))*fxn1*abs(y)/k例4.1.9用计算n阶cotes系数C(n)和求截断误差的主程序，计kb定积分I=af(x)dx的1~3阶Newton- b解（1）先求1~3阶Newton-Cotes的系数和截断误差.输入程>>n1=1,[ 运行后屏幕显示1~3 的系 1, 和截断误差 3n11Cn1 1n2=2Cn2 2n3=3Cn3 3(三)计算Newton-Cotes的程调用格式二调用格式三quad8(‘fun’,a,b,/例4.1.10用梯 计算定积分 解（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')trapzcumtrapzzt、ztcztztcColumns1through Columns250through251 用求积计算定积分.输入程>>symsL=inline(' TS]=quad(L,0,pi/2,1.e-QS FCNTS 用Newton-Cotes求积计算定积分.在6.5中输入程>>symsL=inline(' T8]=quad8(L,0,pi/2,1.e-运行后屏幕显示用Newton-Cotes求积计算定积分的值Q8和递归次数FCNTS分Q8 FCNT8 >>symswzt=abs(Fj-zt),wQS=abs(Fj-QS),wQ8=abs(Fj-Warning:Explicitintegralcouldnotbe>InC:\6p5p1\toolbox\symbolic\@sym\int.matFint(exp(sin(x)),x=0..1/2*pi)Fj= wzt= wQS= wQ8= >>symsxholdonL=inline('Q=quad(L,0,pi/2,1.e-4,2),plot(Q,'g*')holdoff,holdon,h=pi/40;x=0:h:pi/2;y=exp(sin(x));plot(x,ztc,'mp'),holdholdon,[ T8quad8(L,0pi/2,1.e-4,3)holdoffgrid,xlabel('自变量X'),ylabel('因变量Y') 和梯 legend('精确值', 计算定积分','梯形计算定积分4.2（Romberg）及其程4.2.1积分 数值积分 主程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);k=k+1;h=h/2;s=0;forj=1:nx=a+h*(2*j-1);s=s+feval(fun,x);RT(k+1,1)=RT(k,1)/2+h*s;n=2*n;fori=1:k4.2.1取精度为108

Tj1,j1Tj1,

Tj1,j1Tj,

romberg.m

I151dx 01

Tj1,j1Tj1,

（1）取精度分别为108

Tj1,j1Tj1,

symsxfi=int(1/(1+x),x,0,1.5);Fs=double(fi),wR=double(abs(fi-R)),wR1=wR-wugu，wR如下，RTColumns1through0000.95357142857143000.9259835752482800.918741799095710.916297223682920.916905341657170.916290876875400.916444501306570.91629073443494Columns5through000000000 R wugu h Fs= wR= wR1= RT

Tj1,j1Tj,

Columns1through000000 Columns5through 0 0 0 0 0 0 R wuguhwRwR1= RTColumns1through000

Tj1,j1Tj1,0

R wugu h Fs wR wR1

Tj1,j1Tj1,

4.3（Gauss）型积分及其程(一)两点-勒让德积例4.3.1分别用两点–勒让德积分、步长为2的梯形和步长为1的 求积计1解f(x)15

I15xdx，并将计算结果与精确值进行比较.>>x1=-1/sqrt(3);x2=1/sqrt(3);y1=1/(5+x1);y2=1/(5+x2);G2=y1+y2,t=-1:2:1;y=1./(5+t);T=trapz(t,y)symsxFs=double(fi),wQ=double(abs(fi-Q)),wG2=double(abs(fi-G2)),wT=double(abs(fi-T))运行后屏幕显示分别用两点-勒让德积分和步长为2的梯形求积计算I的FswQ依次如下G2 T Fs wQ= wG2= wT=9.044744739561694e-0055.970270275897657e-005(二)n个节点的-勒让德积分及其误差分析的两个程1 计算I1f(x)dx的数值积分及其截断误差 1function[GL,Y,sun=1;su2n=1;su2n1=1;wome=1;symsxfork=1:nwome2=wome^2;Fr=int(wome2,x,-1,1);fork=1:n2syms11RGn=11例4.3.2用–勒让德积分计算I

e2dx3和解>>X1=[-1/(3^(1/2)),1/(3^(1/2))];A1=[1,1];[GL1,Y1,Rn1]=GaussR1(@fun,X1,A1)X2=[-(3/5)^(1/2),0,(3/5)^(1/2)];A2=[5/9,8/9,5/9];[GL2,Y2,Rn2]=GaussR1(@fun,X2,A2)symsFs=double(fi),wGL1=double(abs(fi-GL1)),wGL2=doubleRn1 GL1 Y1 Rn2= GL2= Y2= 451/15261056/2647451/1526Fs= wGL1= wGL2= 1 计算1f(x)dx的数值积分和误差估计的 function[GL,Y,Rn]=GaussR2(fun,X,A,fun2n)n=length(X);n2=2*n;Y=feval(fun,X);1GL=sum(A.*Y);GL=sum(A.*Y);sun=1;su2n=1;su2n1=1;wome=1;symsxforwome2=wome^2;Fr=int(wome2,x,-1,1);fork=1:n2mfun2n=max(fun2n);Rn= 例4.3.3用 计算15xdx，取代数精度为15，再用截断解>>X=[-,,-,- , ,,,,,, ,,x=-1:0.00001:1;fun2n[GL,Y,Rn]=GaussR2(@fun,X,A,fun2n)symstfi=int(1/(5+t),t,-1,1);Fs=double(fi),wGL=double运行后屏幕显示用代数精度为15的-勒让德积分数值计算I的结果GL和精确FswGL依次如下GL Rn0.40546510814876 Y=ColumnsthroughColumnsthroughFswGL在[a,b]上的-勒让德积分及其程b 计算af(x)dx的数值积分 function[GL,Y,Rn]=Gauss(fun,a,b,X,A,fun2n)n=length(X);n2=n*2;T=zeros(1,n);bY=feval(fun,T);GL=((b-a)/2)*sum(A.*Y);sun=1;su2n=1;su2n1=1;wome=1;symsxfork=1:nwome2=wome^2;Fr=int(wome2,x,-1,1);fork=1:n2mfun2n=max(fun2n);Rn=例 利用一般型-勒让德积分计算I

1

e2dx解,,,-,,,,,,,,,,,a=0;b=2;[GL,Y,Rn]=Gausssymsx,fi=int(exp((-x^2)/2)/(sqrt(2*pi)),x,a,b);Fs=double(fi),wGL=double(abs(fi-GL))运行后屏幕显示分别用代数精度为15的一般型-勒让德积分数值计算I的结、、wGLGLRnY Fs —切 求 及 程1x21—切 求 计算I1x21序

f(x)dx及其截断误 主function[ n2=2*n;X=zeros(1,n);A=(pi/n)*ones(1,n);fork=1:nY=feval(fun,X);GC=sum(A.*Y);forsu2n; =(2*pi)*df2n/(211例4.3.5 求 计算I211解>>n=15;[ symsx

fi=int(sqrt(2+x)./sqrt(1-x.^2),x,-1,1);Fs=double(fi)wGL=double(abs(fi-GC))GCXColumns1through