数值积分及其MATLAB实现

第九章 数值积分

9.1 积分的MATLAB符号计算

9.1.1 定积分的MATLAB符号计算

139.1.5 由y?sinx，y?cosx,x??,x?所围成的平面区域D.求平面区域D

22的面积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.36203791318826

9.1.2 变限积分的MATLAB符号计算

例9.1.7 已知F(x)??x2xetsin(2?t3)dt，求F'(x).

解 输入程序：

>> syms x t

F1=int(exp(t)*sin(2+sqrt(t^3)),x,0); F2=int(exp(t)*sin(2+sqrt(t^3)),0,x^2); Fi= F1+ F2; dF=diff(Fi)

运行后屏幕显示计算变限积分F(x)的导数.

9.2 数值积分的思想及其MATLAB程序

9.2.3 矩形公式的MATLAB程序

（一） 函数sum的调用格式

调用格式一：sum(X)

调用格式二：sum (X,DIM)

例9.2.2 用MATLAB和矩形公式（9.3）、（9.4）计算较.

解 将[0,?/2]分成20等份，步长为?/40，输入程序

>> 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

129 ．

??20e

sinxdx，并与精确值比

第九章 数值积分

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）计算

??20e

?xsinxdx，并与精确值比较.

解 将[0,?/2]分成20等份，步长为?/40，输入程序如下（注意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.3961

wz1 = wz2 =

0.0088 0.0075

9.3 插值型数值积分及其MATLAB 程序

9.3.2 梯形公式的MATLAB程序

(一) 根据梯形公式和估计误差公式自己编写MATLAB程序计算定积分

例9.3.2 分别取h??/8000,?/800,?/80，用梯形公式计算定积分

I??响.

?/20e

sinxdx，并与精确值比较.然后观察h对计算结果的有效数字和绝对误差的影

解 编写并输入如下程序

>>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 t

f=exp(sin(t)); intf=int(f,t,a,b), Fs=double(intf), Juewucha8000=abs(z8000-Fs)

运行后屏幕显示取h??/8000时，积分区间[0,?/2]上等距节点的个数n，用梯形公式计算定积分I的值z8000和精确值intf的近似值Fs及其绝对误差Juewucha8000.

(二) 用函数trapz计算定积分

调用格式一：Z =trapz(Y) 调用格式二：Z =trapz(X,Y)

130 ．

第九章 数值积分

调用格式三：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分别计算

?4?

?/20

e

?xsinxdx，精确

到10，并与矩形公式（9.3），（9.4）比较.

解 将[0,?/2]分成20等份，步长为?/40，输入程序如下（注意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)/2 z3=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);

end

T(i+1)=T(i)/2+h*s; end

T=T(1:m);

例9.3.6 用rctrap计算I?1?2??20e

?x22dx，递归14次，并将计算结果与精确

值比较.

解 输入程序

>>T=rctrap(@fun,0,pi/2,14), syms t

fi=int(exp((-t^2)/2)/(sqrt(2*pi)),t,0, pi/2); Fs= double(fi), wT= double(abs(fi-T))

运行后屏幕显示I精确值Fs，用rctrap计算I的递归值T和T与精确值Fs的绝对误差wT

9.3.3 辛普森公式及其误差分析

?edx，取n?20001个等距节点，并

2?0将计算结果与精确值比较，然后再取n?13计算，观察n对误差的影响.

解 由n?2m?1?20001，得m?10000.根据辛普森（Simpson）公式编写并输

例9.3.7 用辛普森公式计算I?入下面的程序

>> 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((-t^2)/2)/(sqrt(2*pi)); intf=int(f,t,a,b), Fs=double(intf); Juewucha=abs(z-Fs)

运行后屏幕显示用辛普森公式(9.11)计算定积分I的近似值z和精确值intf及其绝对误差Juewucha（取n?20001个等距节点）.

131 ．

11?x22第九章 数值积分

例9.3.8 估计用辛普森公式计算定积分I???20e

sinxdx时的误差，取h??/40.

解 根据估计误差公式，先输入求f(4)(x)的程序

>>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=(h^4)*(pi/2)*max(juyx4)/180

运行后屏幕显示误差估计值

RS =

3.610450295892220e-006

9.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); end

for k=3:2:2*m

x3=x+k*h; z3= z3+4*feval (fun,x3); end

y=(z1+z2+z3)*h/3; 例9.3.9 用comsimpson.m和quad.m分别计算定积分I??41?2?10e

?x22dx，取

精度为10，并与精确值比较.

解 输入程序

>> [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) )

运行后屏幕显示I的精确值Fs，用comsimpson.m和quad.m分别计算I的近似值Q2、Q1和迭代次数FCNT14，取精度分别为10，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.0745230082 Q1 = FCNT14 = Q2 =

0.3413 13 0.3413 Fs = wQ1 = wQ2 =

0.3413 3.6619e-009 3.7061e-005

?4

132 ．

第九章 数值积分

9.3.6 牛顿-科茨（Newton-Cotes）数值积分和误差分析的MATLAB程序

(一) 估计误差的MATLAB程序

计算n阶牛顿-科茨的公式的截断误差公式的MATLAB主程序

function RNC=ncE(n) suk=1; p=n/2-fix(n/2); if p==0

for k=1:n+2 suk=suk*k; end

suk; syms t a b fxn2,su=t^2; for u=1:n

su=su*(t-u); end

su; intf=int(su,t,0,n); y=double(intf); RNC= (((b-a)/n)^(n+3))*fxn2*abs(y)/ suk; else

for k=1:n+1 suk=suk*k; end

suk; syms t a b fxn1,su=t; for u=1:n

su=su*(t-u); end

su; intf=int(su,t,0,n); y=double(intf); RNC= (((b-a)/n)^(n+2))*fxn1*abs(y)/ suk; end

例9.3.14 用求截断误差公式的MATLAB主程序，求计算定积分值的1,2,3,4,8阶牛顿–科茨公式的截断误差公式.

解 输入求1,2,3,4,8阶牛顿–科茨公式的截断误差公式的程序

>> 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*fxn1 n = RNC2 =

2 1/90*(1/2*b-1/2*a)^5*fxn2 n = RNC3 =

3 3/80*(1/3*b-1/3*a)^5*fxn1 n = RNC4 =

4 8/945*(1/4*b-1/4*a)^7*fxn2 n = RNC8 =

8 35065906189543/6926923254988800*(1/8*b-1/8*a)^11*fxn2

例9.3.15 用MATLAB软件估计用5、6阶牛顿–科茨公式计算定积分I??baf(x)dx的近似

??20e

sinxdx的截断误差公式和误差限，取15位小数.

解 输入求n?5,6阶牛顿–科茨公式的截断误差公式和被积函数的6，8阶导函数的程序

>> n=5;RNC5=ncE(n), n=6;RNC6=ncE(n) syms x

y=exp(sin(x)); yx6=diff(y,x,6), yx8=diff(y,x,8)

运行后输出被积函数的6，8阶导函数（略）和n?5,6阶牛顿–科茨公式的截断误差公式

133 ．

第九章 数值积分

如下

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*myx6

RNC6 =9/1400*(1/6*b-1/6*a)^9*myx8

运行后屏幕显示误差限如下

RNC5 = RNC6 =

3.625385713339320e-004 3.826183594914823e-005

(n)(二) 计算科茨系数Ck和求截断误差公式的MATLAB程序 (n)计算n阶科茨系数Ck和求截断误差公式的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:n

su1=su*(t-u);m01=m0*u; su=su1;m0=m01; end

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; for j=1:n

k1=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:divideByZero end

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==0

for k=1:n+2 suk=suk*k; end

suk; syms t a b fxn2,su=t^2; for u=1:n

su=su*(t-u); end

su; intf=int(su,t,0,n); y=double(intf); RNCn=(((b-a)/n)^(n+3))*fxn2*abs(y)/suk; else

for k=1:n+1 suk=suk*k;

134 ．

第九章 数值积分

end

suk; syms t a b fxn1,su=t; for u=1:n

su=su*(t-u); end

su; intf=int(su,t,0,n); y=double(intf); RNCn=(((b-a)/n)^(n+2))*fxn1*abs(y)/ suk; end

(n)例9.3.16 用计算n阶科茨系数Ck和求截断误差公式的MATLAB主程序，计算

定积分I=

?baf(x)dx的1~3阶牛顿–科茨公式的系数和截断误差公式.

解 （1）先求1~3阶牛顿–科茨公式的系数和截断误差公式.输入程序

>> n1=1,[Cn1,RNCn1]=newcotE(n1) n2=2,[Cn2,RNCn2]=newcotE(n2) n3=3,[Cn3,RNCn3]=newcotE(n3)

运行后屏幕显示1~3阶牛顿–科茨公式的系数Cn1, Cn2, Cn3和截断误差公式RNCn1, RNCn2, RNCn3如下

n1 =

1 Cn1 =

1/2 1/2 RNCn1 =

1/12*(b-a)^3*fxn1 n2 =

2 Cn2 =

1/6 2/3 1/6 RNCn2 =

1/90*(1/2*b-1/2*a)^5*fxn2 n3 =

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 用梯形公式、辛普森和牛顿–科茨求积公式计算定积分

I???/20e

sinxdx，取精度10，作出它们的积分图，并与精确值进行比较;

?4

解 （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.10437572798742 ztc =

Columns 1 through 3

135 ．

第九章 数值积分

0 0.00630298652792 0.01264569951380 ……………………………………………………………………….. Columns 250 through 251

3.08729642810745 3.10437572798742 （2）用辛普森求积公式计算定积分. 输入程序

>> syms x

L= inline(' exp(sin(x))');

[QS,FCNTS] =quad(L,0, pi/2,1.e-4,2)

运行后屏幕显示用辛普森求积公式计算定积分的值QS和递归次数FCNTS分别如下

QS = FCNTS =

3.10438133817254 13

（3）用牛顿–科茨求积公式计算定积分. 在MATLAB6.5中输入程序

>> syms x

L= inline(' exp(sin(x))');

[Q8,FCNT8] = quad8(L,0, pi/2,1.e-4,3)

运行后屏幕显示用牛顿–科茨求积公式计算定积分的值Q8和递归次数FCNTS分别如下

Q8 = FCNT8 = 3.10437901785555 33 （4）输入求定积分的精确值的程序

>> syms x

y=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:\\MATLAB6p5p1\\toolbox\\symbolic\\@sym\\int.m at line 58 F =

int(exp(sin(x)),x = 0 .. 1/2*pi) Fj = wzt =

3.10437901785556 3.289868133471430e-006 wQS = wQ8 =

2.320316987436399e-006 7.993605777301127e-015 （5）输入画图程序

>> syms x

F=int(exp(sin(x)),0,pi/2),Fj=double(F),plot(Fj,'ro'), hold on

L= 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 off

hold on, [Q8,FCNT8] = quad8(L,0, pi/2,1.e-4,3), hold off grid, xlabel('自变量 X'), ylabel('因变量 Y') title('牛顿–科茨公式和梯形公式计算数值积分') legend('精确值', ' 辛普森公式计算定积分','梯形公式计算定积分','牛

顿–科茨公式计算定积分') 运行后屏幕显示图形（略）.

9.3.7 利用三次样条求表格型数值积分的MATLAB方法

1?x2例9.3.18 给出概率积分f(x)?e 的数据表： 2?X f(X)

136 ．

0.46 0.47 0.48 0.49 0.322 9 0.319 9 0.316 8 0.313 8 第九章 数值积分

试首先用三次样条构造被积函数P(x)，再分别用quad和quadl函数计算定积分

I??0.490.46f(t)dt，精度为10?4，并将计算结果与精确值比较.

解 方法1 首先用三次样条构造被积函数P(x)，再分别用quad和quadl函数计算定积分I.在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)

运行后屏幕显示分别用辛普森和牛顿–科茨公式计算I的结果为

pp =

form: 'pp'

breaks: [0.4600 0.4700 0.4800 0.4900] coefs: [3x4 double] pieces: 3 order: 4 dim: 1

Q = FCNT = Ql = FCNTl =

0.0096 13 0.0096 18 方法2 计算I的精确解y及其近似解y1. 输入程序

>> syms x

f=exp(-x.^2)./( sqrt(2*pi)); F= int(f, 0.46, 0.49) y=simple(F), y1=double(y)

运行后屏幕显示结果

y =

1125899906842624/5644425081792261*pi^(1/2)*(erf(49/100)-erf(23/50)) y1 =

0.0096

方法3 首先用三次样条构造被积函数P(x)，再分别用辛普森和牛顿–科茨公式计算定积分I.在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 13

9.3.8 利用拉格朗日插值等方法求表格型数值积分的MATLAB方法

例9.3.19 测得定积分I???20f(x)dx中被积函数Y?f(x)在积分变量x的某几

个特定点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函数计算，取精度为10，并将计算结果与精确值1进行比较.

解 方法1 (1) 拉格朗日插值多项式的MATLAB主程序

function L=lfun(X,Y)

137 ．

?4第九章 数值积分

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)); end end

l(k,:)=poly2sym (V); end

L=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 x

L1=inline('-14644281526214207/140737488355328000*x.^2+22

80009553133670687/2251799813685248000*x');

L2=inline('-88810055001957943/351843720888320000*x.^2+15

892292661979563/14073748835532800*x-15511422992738729/703687441776640000');

L3=inline('-261101353037957127/703687441776640000*x.^2+3

62054074498830351/281474976710656000*x-207985426949036897/2814749767106560000');

L4=inline('-131688789655946847/281474976710656000*x.^2+1

04193434945799191/70368744177664000*x-9309740592764125768759/54295819320043765760000');

[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+ FCNT44

[Q841,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函数计算I的近似值Q、Q8和迭代次数FCNT、FCNT8，取精度为10如下

Q = FCNT = Q8 = FCNT8 =

0.99957595925413 52 0.99957595925413 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 x

138 ．

?4

第九章 数值积分

L1=inline('-9070382563990323/56294995342131200*x.^3-793956227651389/281474976710656000*x.^2+2253151961697736691/2251799813685248000*x');

L2=inline('-41118170146065463/351843720888320000*x.^3-17

00780706742359/21990232555520000*x.^2+735208215754549839/703687441776640000*x-25688208496779769/2814749767106560000');

L3 =

inline('-131688789655946847/281474976710656000*x.^2+1041

93434945799191/70368744177664000*x-9309740592764125768759/54295819320043765760000');

[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+ FCNT43

[Q841,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函数计算I的近似值Q、Q8和迭代次数FCNT、FCNT8，取精度为10如下

Q = FCNT=

0.99975250938584 39 Q84= FCNT8=

0.99975250938584 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.00012965800727 33 1.00012956660843 17

?49.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((wugu>wucha)&(k

k=k+1; h=h/2; s=0; for j=1:n

x=a+h*(2*j-1); s=s+feval(fun,x);

end

RT(k+1,1)= RT(k,1)/2+h*s; n=2*n; for i=1:k

RT(k+1,i+1)=((4^i)*RT(k+1,i)-RT(k,i))/(4^i-1); end

wugu=abs(RT(k+1,k)-RT(k+1,k+1)); end

139 ．

第九章 数值积分

R=RT(k+1,k+1);

?8例9.4.2 取精度为??10，分别用Tj?1,j?1?Tj?1,j??和Tj?1,j?1?Tj,j??作

1?01?xdx，龙贝格积分表，最小步

?8?7长h，取精度为10,并与精确值比较.然后取精度为10，用Tj?1,j?1?Tj?1,j??作为计

为计算停止的条件，用romberg.m程序计算I?1.5算停止的条件，观察用龙贝格求积公式计算的结果与精确值的绝对误差是否满足10.

?8解 （1）取精度分别为10，当用Tj?1,j?1?Tj?1,j??作为计算停止的条件，输入

?7程序

>>F=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

?8运行后屏幕显示I精确值Fs，取精度为10，利用龙贝格求积公式的romberg.m程序计算I的近似值R，龙贝格积分表RT，误差估计err，最小步长h，近似值R与精确值Fs的绝对误差wR如下

RT =

Columns 1 through 4

1.05000000000000 0 0 0 0.95357142857143 0.92142857142857 0 0 0.92598357524828 0.91678762414057 0.91647822765470 0 0.91874179909571 0.91632787371152 0.91629722368292 0.9162943506040 0.91690534165717 0.91629318917766 0.91629087687540 0.91629077613242 0.91644450130657 0.91629088785636 0.91629073443494 0.91629073217398

Columns 5 through 6

0 0 0 0 0 0 0 0

0.91629076211489 0 0.91629073200160 0.91629073197216

R = wugu = h =

0.91629073197216 2.943623123030648e-011 0.04687500000000 Fs = wR = wR1 =

0.91629073187416 9.800686628988125e-011 6.857063505957478e-011

如果用Tj?1,j?1?Tj,j??作为计算停止的条件，运行的结果如下

RT =

Columns 1 through 4

1.05000000000000 0 0 0 0.95357142857143 0.92142857142857 0 0 0.92598357524828 0.91678762414057 0.91647822765470 0

0.91874179909571 0.91632787371152 0.91629722368292 0.91629435060401 0.91690534165717 0.91629318917766 0.91629087687540 0.91629077613242 0.91644450130657 0.91629088785636 0.91629073443494 0.91629073217398 0.91632918157305 0.91629074166188 0.91629073191558 0.91629073187559

Columns 5 through 7

0 0 0 0 0 0 0 0 0 0 0 0 0.91629076211489 0 0 0.91629073200160 0.91629073197216 0 0.91629073187442 0.91629073187429 0.91629073187427 R = wugu = h =

0.91629073187427 9.789258292869363e-011 0.02343750000000 wR = wR1 =

1.142833611876187e-013 -9.777829956750600e-011

140 ．

第九章 数值积分

?7（2）取精度分别为10，当用Tj?1,j?1?Tj?1,j??作为计算停止的条件，计算结果

如下

RT =

Columns 1 through 4

1.05000000000000 0 0 0

0 0 0 0.91629076211489

R = wugu = h =

0.91629076211489 1.401753557672691e-008 0.09375000000000 Fs = wR = wR1 =

0.91629073187416 3.024073362534249e-008 1.622319804861558e-008

?7由此可见，取精度为10，用Tj?1,j?1?Tj?1,j??作为计算停止的条件，用龙贝格

求积公式计算的结果与精确值的绝对误差满足10.

?7

9.5 自适应积分及其MATLAB程序

例9.5.1 用quad.m程序计算

?1.501?7dx，取精度为10, 并与精确值作比较. 1?x解 输入程序

>> F = inline('1./(1+x)'); [Q,FCNT] = quad(F,0,1.5,1.e-7,3) syms x

fi=int(1/(1+x),x,0,1.5);

Fs=double (fi), wQ= double (abs(fi-Q))

?7运行后屏幕显示I精确值Fs，取精度为10，用函数quad计算I的近似值Q，递归次数FCNT它与精确值Fs的绝对误差wQ如下

Q = FCNT =

0.91629073716847 29 Fs = wQ =

0.91629073187416 5.294317803384376e-009

9.6 高斯（Gauss）型积分公式及其MATLAB程序

9.6.2 在[?1,1]上的高斯-勒让德积分及其MATLAB程序

(一) 两点高斯-勒让德积分公式

例9.6.3 分别用两点高斯–勒让德积分公式（9.60）、步长为2的梯形公式和步长为1的辛普森求积公式计算I?解 设f(x)?1??15?xdx，并将计算结果与精确值进行比较.

11. 输入程序 5?x>> 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) syms x

fi=int(1/(5+x),x,-1,1);

Fs= double (fi), Q=0.40555555555556;

wQ= double (abs(fi-Q)), wG2= double (abs(fi-G2)), wT= double (abs(fi-T))

运行后屏幕显示分别用两点高斯-勒让德积分公式和步长为2的梯形求积公式计算I的结果G2、T 、Q和精确值Fs及其分别与Fs的绝对误差w G2、wT和步长为1的辛普森与Fs的绝对误wQ依次如下

G2 = T = Fs =

141 ．

第九章 数值积分

0.40540540540541 0.41666666666667 0.40546510810816 wQ = wG2 = wT =

9.044744739561694e-005 5.970270275897657e-005 0.01120155855850

(二) n个节点的高斯-勒让德积分公式及其误差分析的两个MATLAB程序 用高斯-勒让德积分公式计算I?主程序

function [GL,Y, RGn]=GaussR1(fun,X,A)

n=length(X);n2=2*n;Y=feval(fun,X);GL=sum(A.*Y); sun=1; su2n=1; su2n1=1; wome=1; syms x for k=1:n

wome=wome*(x-X(k)); end

wome2= wome^2;Fr=int(wome2,x,-1,1); for k=1:n2

su2n=su2n*k; end syms M

RGn= Fr*M/su2n;

例9.6.4 用高斯–勒让德积分公式计算I??1?1f(x)dx的数值积分及其截断误差公式的MATLAB

1?2?1?1e

?x22dx，取代数精度为3和5，

再根据截断误差公式（9.62）写出误差公式，并将计算结果与精确值进行比较.

解 输入程序

>> 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) syms x

fi=int(exp( (-x.^2)./2)./(sqrt(2*pi)),x,-1,1);

Fs= double (fi), wGL1 = double (abs(fi-GL1)), wGL2= double

(abs(fi-GL2)) 运行后屏幕显示如下

Rn1 = GL1 = Y1 =

1/135*M 1754/2597 877/2597 877/2597 Rn2 = GL2 = Y2 =

1/15750*M 237/347 451/1526 1056/2647 451/1526 Fs = wGL1 = wGL2 =

2315/3391 107/14668 36/116971

用高斯-勒让德积分公式计算

?1?1f(x)dx的数值积分和误差估计的MATLAB主程序

function [GL,Y,Rn]=GaussR2 (fun,X,A,fun2n) n=length(X);n2=2*n; Y=feval(fun,X); GL=sum(A.*Y); GL=sum(A.*Y); sun=1; su2n=1; su2n1=1; wome=1; syms x for k=1:n

wome=wome*(x-X(k)); end

wome2= wome^2;Fr=int(wome2,x,-1,1); for k=1:n2

su2n=su2n*k; end

mfun2n =max(fun2n); Rn = Fr*mfun2n/su2n;

142 ．

第九章 数值积分

例9.6.5 用高斯–勒让德积分公式计算

1??15?xdx，取代数精度为15，再用截断

1误差公式（9.62）估计误差，并将计算结果与精确值进行比较. 解 输入程序

>> X=[-0.9602898565, -0.7966664774, -0.5255324099,

-0.1834346425, 0.9602898565, 0.7966664774, 0.5255324099, 0.1834346425];

A=[0.1012285363, 0.2223810345, 0.3137066459, 0.3626837834,

0.1012285363, 0.2223810345, 0.3137066459, 0.3626837834];

x=-1:0.00001:1; fun2n =20922789888000./(5+x).^17; [GL,Y,Rn]=GaussR2 (@fun,X,A,fun2n) syms t

fi=int(1/(5+t),t,-1,1); Fs= double (fi), wGL = double

(abs(fi-GL))

运行后屏幕显示用代数精度为15的高斯-勒让德积分公式数值计算I的结果GL和精确值Fs及其两者的绝对误差wGL依次如下

GL = Rn =

0.4 0546510814876 2.709468201849595e-015 Y =

Columns 1 through 4

0.24754251282336 0.23790641276104 0.22349027674545 0.20761682356139

Columns 5 through 8

0.16777707528929 0.17251294410310 0.18097803538503 0.19292227431611

Fs = wGL =

0.40546510810816 4.059186608518054e-011

9.6.3 在[a,b]上的高斯-勒让德积分公式及其MATLAB程序

用一般型高斯-勒让德积分公式计算

?baf(x)dx的数值积分MATLAB主程序

function [GL,Y,Rn]=Gauss (fun,a,b,X,A,fun2n) n=length(X); n2=n*2; T=zeros(1,n); T=((a+b)/2)+((b-a)/2)*X;

Y=feval(fun,T); GL=((b-a)/2)*sum(A.*Y); sun=1; su2n=1; su2n1=1; wome=1; syms x for k=1:n

wome=wome*(x-X(k)); end

wome2= wome^2;Fr=int(wome2,x,-1,1); for k=1:n2 su2n=su2n*k; end

mfun2n =max(fun2n); Rn = Fr*mfun2n/su2n;

例9.6.6 利用一般型高斯-勒让德积分公式计算I?1?2?20e

?x22dx，取代数精度

为15，并与精确值比较.

解 输入程序

>>X=[-0.9602898565, -0.7966664774, -0.5255324099, -0.1834346425, 0.9602898565, 0.7966664774, 0.5255324099, 0.1834346425];

A=[0.1012285363, 0.2223810345, 0.3137066459, 0.3626837834, 0.1012285363, 0.2223810345, 0.3137066459, 0.3626837834];

143 ．

第九章 数值积分

a=0;b=2; [GL,Y,Rn]=Gauss (@fun,a,b,X,A,fun2n) syms x, fi=int(exp((-x^2)/2)/(sqrt(2*pi)),x,a,b); Fs= double (fi), wGL = double (abs(fi-GL))

运行后屏幕显示分别用代数精度为15的一般型高斯-勒让德积分公式数值计算I的结果GL、高斯-勒让德积分公式的误差Rn、估计被积函数值Y和精确值Fs及其两者的绝对误差wGL依次如下

GL = Rn =

0.47724986810026 2.709468201849595e-015 Y =

0.39862785922733 0.39077989968078 0.35647246460968 0.28583858238797 0.05840774991356 0.07942486895803 0.12461041358066 0.19805761421327

Fs = wGL =

0.47724986805182 4.843847992335364e-011

9.6.4 拉道积分公式和洛巴托积分公式及其MATLAB程序

(一) 用洛巴托数值积分的MATLAB程序

调用格式一：Q = quadl (FUN,a,b)

调用格式二：Q = quadl (FUN,a,b,TOL) 调用格式三：[Q,FCNT] = quadl (...)

调用格式四：quadl (FUN,a,b,TOL,TRACE)

调用格式四：quadl (FUN,A,B,TOL,TRACE,P1,P2,...) 调用格式五：quadl FUN(X,P1,P2,...)

例9.6.7 用函数quadl和quad8及quad分别数值计算I??6?141?2??20e

?x22dx，分

别取精度为10和10，将计算结果分别与精确值比较，并画出图形.

解 输入程序

>> subplot(2,1,1),

[Ql,FCNTl]=quadl(@fun,0,pi/2,1.e-6) [Q8,FCNT8]=quad8(@fun,0,pi/2,1.e-6,3) [Q,FCNT]=quad(@fun,0,pi/2,1.e-6) syms x

fi=int(exp((-x.^2)./2)./(sqrt(2*pi)),x,0, pi/2); Fs=double(fi),

wQl=double(abs(fi-Ql)) wQ8=double(abs(fi-Q8)), wQ=double(abs(fi-Q))

xlabel('自变量 X'), ylabel('因变量 Y')

title('取精度为1.e-6，被积函数f(x)=exp(-x^2/2)/(sqrt(2 pi))在[0,pi/2]上的数值积分')

legend('取精度1.e-6，用牛顿–科茨公式计算积分过程的散点图') hold on

subplot(2,1,2)

[Ql14,FCNTl14]=quadl(@fun,0,pi/2,1.e-14) [Q814,FCNT814]=quad8(@fun,0,pi/2,1.e-14,3) [Q14,FCNT14]=quad(@fun,0,pi/2,1.e-14,1) wQl14=double(abs(fi-Ql14)), wQ814=double(abs(fi-Q814)) wQ14=double(abs(fi-Q14)),

xlabel('自变量 X'), ylabel('因变量 Y')

title('取精度为1.e-14，被积函数f(x)=exp(-x^2/2)/(sqrt(2 pi))在[0,pi/2]上的数值积分')

legend('取精度1.e-14，用牛顿–科茨公式计算积分过程的散点图') hold off

144 ．

第九章 数值积分

运行后屏幕显示如下

Ql = FCNTl = Q8 = FCNT8 = 0.44188501970377 18 0.44188501721659 33 Q = FCNT = Fs =

0.44188511819589 13 0.44188501721659 wQl = wQ8 =

2.487183090649321e-009 3.572753311048285e-015 wQ = Ql14 =

1.009793004969618e-007 0.44188501721659, FCNTl14 = Q814 = FCNT814 = 138 0.44188501721659 129 Q14 = FCNT14 = wQl14 =

0.44188501721659 553 1.310619347103004e-016 wQ814 = wQ14 =

9.098267021473094e-017 2.003963224778471e-017

例9.6.8 用辛普森求积公式、牛顿–科茨求积公式和洛巴托求积公式分别计算

I???20excos2xdx，取精度分别为103?14和10，并与精确值比较.请再改用取10个

?4?4节点，分别用三阶拉格朗日插值多项式和三次样条构造被积函数，然后用辛普森求积公式、牛顿–科茨求积公式和洛巴托求积公式分别计算I，取计算精度分别为10.再改用100, 1 000,10 000个节点，用三次样条构造被积函数计算，观察计算结果与精确值的绝对误差有何变化.

解 （1）建立名为fun.m的M文件

function y = fun(x)

y = exp(x.^3).*cos(2*x);

?4取精度为10，输入程序

>> F=inline('exp(x.^3).*cos(2*x)'); [Ql,FCNTl]=quadl(@fun,0,pi/2,1.e-4,3) [Q8,FCNT8]=quad8(@fun,0,pi/2,1.e-4,3) [Q,FCNT]=quad(@fun,0,pi/2,1.e-4,3) syms x

fi=int(exp(x.^3).*cos(2*x),x,0, pi/2); Fs=double(fi), wQl=double(abs(fi-Ql) )

wQ8=double(abs(fi-Q8) ), wQ=double(abs(fi-Q) )

运行后屏幕显示如下

Ql =

-6.45505792320975 FCNTl = 48 Q8 =

-6.45505883179317 FCNT8 = 33 Q =

-6.45506213501363 FCNT = 33

Warning: Explicit integral could not be found. > In C:\\MATLAB6p5\\toolbox\\symbolic\\@sym\\int.m at line 58 Fs =

-6.45505792116786 wQl =

2.041896319084332e-009 wQ8 =

9.106253137107240e-007 wQ =

145 ．

第九章 数值积分

4.213845774426019e-006 （2）取精度为10，输入程序

>> F=inline('exp(x.^3).*cos(2*x)'); [Ql,FCNTl]=quadl(@fun,0,pi/2,1.e-14,3) [Q8,FCNT8]=quad8(@fun,0,pi/2,1.e-14,3) [Q,FCNT]=quad(@fun,0,pi/2,1.e-14,3) syms x

fi=int(exp(x.^3).*cos(2*x),x,0, pi/2); Fs=double(fi), wQl=double(abs(fi-Ql))

wQ8=double(abs(fi-Q8)), wQ=double(abs(fi-Q))

运行后屏幕显示如下

Ql =

-6.45505792116785 FCNTl = 798 Q8 =

-6.45505792116785 FCNT8 =

4273 Q =

-6.45505792116785 FCNT =

2961

Warning: Explicit integral could not be found. > In C:\\MATLAB6p5\\toolbox\\symbolic\\@sym\\int.m at line 58 Fs =

-6.45505792116786 wQl =

3.663345763769995e-015 wQ8 =

1.886988924369744e-015 wQ =

2.775167344069869e-015

这三种方法计算结果的绝对误差都小于10，达到所要求的精度.

（3）取10个节点，改用三阶拉格朗日插值多项式构造被积函数，然后用辛普森求积公式、牛顿–科茨求积公式和洛巴托求积公式分别计算I，取精度分别为10，输入程序

>>X=[0,0.1745,0.3491,0.5236,0.6981,0.8727,1.0472,1.2217,

1.3963,1.5708];

Y= exp(X.^3).*cos(2* X); L1=lfun(X(1:4),Y(1:4)), L2=lfun(X(4:7),Y(4:7)), L3=lfun(X(7:10),Y(7:10)),

运行后输出多项式L1，L2，L3.然后输入程序

>> syms x L1 = inline('

528475254238581253770784020141/1267650600228229401496703205376*x.^3-2150921733462980567549357641495/1267650600228229401496703205376*x.^2-110805035875068143559156492912389/3318709271397504573118368991674368*x+1');

L2 = inline('

-32514595309399085335529842115599/2535301200456458802993406410752*x.^3+28880777683589964816385504370251/1267650600228229401496703205376*x.^2-9690424145215092365005571145155/633825300114114700748351602688*x+21181277275320971857213583641339/5070602400912917605986812821504');

L3 = inline('

-22347657567307906685985592460479/39614081257132168796771975168*x.^3+19439764627825105995989937383455/9903520314283042199192993792*x.^2-45346855519117547561405281365707/19807040628566084398

146 ．

?14?14?14第九章 数值积分

385987584*x+17651524197374204698139309247929/19807040628566084398385987584');

[Q1,FCNT1] = quad(L1,0, 0.5236,1.e-4);

[Q2,FCNT2] =quad(L2, 0.5236, 1.0472,1.e-4); [Q3,FCNT3] = quad(L3, 1.0472, pi/2,1.e-4); Q= Q1+ Q2+ Q3 ,FCNT= FCNT1+ FCNT2+ FCNT3 [Q81,FCNT81] = quad8(L1,0, 0.5236,1.e-4);

[Q82,FCNT82] =quad8(L2, 0.5236, 1.0472,1.e-4); [Q83,FCNT83] = quad8(L3, 1.0472, pi/2,1.e-4); Q8= Q81+ Q82+ Q83 ,FCN8T= FCNT81+ FCNT82+ FCNT83 [Ql1,FCNTl1] = quadl(L1,0, 0.5236,1.e-4);

[Ql2,FCNTl2] =quadl(L2, 0.5236, 1.0472,1.e-4); [Ql3,FCNTl3] = quadl(L3, 1.0472, pi/2,1.e-4); Ql= Ql1+ Ql2+ Ql3 ,FCNTl= FCNTl1+ FCNTl2+ FCNTl3 syms x

fi=int(exp(x.^3).*cos(2*x),x,0, pi/2);

Fs= double (fi), wQl= double (abs(fi-Ql) )

wQ8= double (abs(fi-Q8) ), wQ= double (abs(fi-Q) )

运行后屏幕显示如下

Q = FCNT= Q8 = FCN8T =

-6.63556515894584 39 -6.63556515894593 99 Ql = FCNTl = -6.63556515894607 54

Warning: Explicit integral could not be found.

> In D:\\MATLAB6P5\\toolbox\\symbolic\\@sym\\int.m at line 58 Fs = wQl =

-6.45505792116786 0.18050723777821 wQ8 = wQ =

0.18050723777807 0.18050723777798

（4）取10个节点，改用三次样条构造被积函数，然后用辛普森求积公式、牛顿–科茨求积公式和洛巴托求积公式分别计算I，取精度分别为10.

首先建立并保存名为scjfQl8.m的M文件

function wQ=scjfQl8(n) a=0; b= pi/2;X=0:pi/(2*n-1):pi/2; fs= exp(X.^3).*cos(2* X); pp=spline(X,fs); Qsl=quadl(@ppval,a,b,[],[],pp) Qs8=quad8(@ppval,a,b,[],[],pp), Qs=quad(@ppval,a,b,[],[],pp) syms x

fi=int(exp(x.^3).*cos(2*x),x,0, pi/2); Fs=double(fi) wQl=double(abs(fi-Qsl)), wQ8=double(abs(fi-Qs8)) wQ=double(abs(fi-Qs)) 然后输入程序

>> n=10;wQ=scjfQl8(n)

运行后屏幕显示结果如下

Qsl = Qs8 = Qs = -6.33658113206801 -6.33658632890273

-6.33659053214403

Warning: Explicit integral could not be found. > In C:\\MATLAB6p5\\toolbox\\symbolic\\@sym\\int.m at line 58 In C:\\MATLAB6p5\\work\\scjfQl8.m at line 7 Fs =

-6.45505792116786

wQl = wQ8 = wQ =

0.11847678909985 0.11847159226513 0.11846738902383

?4 147 ．

第九章 数值积分

例9.6.10 测得定积分I???20f(x)dx中被积函数Y?f(x)在积分变量x的某几

个特定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，试根据这些值首先分别用分段二阶、三阶拉格朗日插值多项式和三次样条插值构造被积函数，然后分别用辛普森和洛巴托求积公式计算，取精度分别为10?14和

10?4，并与精确值1进行比较.

解 下面分别用分段二阶、三阶拉格朗日插值多项式和三次样条插值方法分别构造被积函数，计算数值积分.

方法1 (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: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，然后输入程序

>> syms x

L1=inline('-14644281526214207/140737488355328000*x.^2+22

80009553133670687/2251799813685248000*x');

L2=inline('-88810055001957943/351843720888320000*x.^2+15

892292661979563/14073748835532800*x-15511422992738729/703687441776640000');

L3=inline('-261101353037957127/703687441776640000*x.^2+3

62054074498830351/281474976710656000*x-207985426949036897/2814749767106560000');

L4=inline('-131688789655946847/281474976710656000*x.^2+1

04193434945799191/70368744177664000*x-9309740592764125768759/54295819320043765760000');

[Q141,FCNT141] = quad(L1,0, 0.4712,1.e-14);

[Q142,FCNT142] =quad(L2, 0.4712, 0.6283,1.e-14); [Q143,FCNT143] = quad(L3, 0.6283, 1.0210,1.e-14); [Q144,FCNT144] = quad(L4, 1.0210, pi/2,1.e-14); Q14= Q141+ Q142+ Q143+ Q144,

FCNT14= FCNT141+ FCNT142+ FCNT143+ FCNT144 [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); Q4= Q41+ Q42+ Q43+ Q44,

FCNT4= FCNT41+ FCNT42+ FCNT43+ FCNT44

[Ql41,FCNTl41] = quadl(L1,0, 0.4712,1.e-4);

[Ql42,FCNTl42] =quadl(L2, 0.4712, 0.6283,1.e-4); [Ql43,FCNTl43] = quadl(L3, 0.6283, 1.0210,1.e-4); [Ql44,FCNTl44] = quadl(L4, 1.0210, pi/2,1.e-4); Ql4= Ql41+ Ql42+ Ql43+ Ql44,

FCNTl4= FCNTl41+ FCNTl42+ FCNTl43+ FCNTl44

[Ql141,FCNTl141] = quadl(L1,0, 0.4712,1.e-14,2);

[Ql142,FCNTl142] =quadl(L2, 0.4712, 0.6283,1.e-14,2); [Ql143,FCNTl143] = quadl(L3, 0.6283, 1.0210,1.e-14,2); [Ql144,FCNTl144] = quadl(L4, 1.0210, pi/2,1.e-14,2); Ql14= Ql141+ Ql142+ Ql143+ Ql144,

FCNTl14= FCNTl141+ FCNTl142+ FCNTl143+ FCNTl144 wQ14=double(abs(1- Q14)), wQ4=double(abs(1- Q4))

wQl4=double(abs(1- Ql4)), wQl14=double(abs(1- Ql14))

运行后屏幕显示如下

Q14 = FCNT14 = Q4 = FCNT4 =

148 ．

第九章 数值积分

0.99957595925413 52 0.99957595925413 52 Ql4 = FCNTl4 = Ql14 = FCNTl14 = 0.99957595925413 72 0.99957595925413 72 wQ14 = wQ4 =

4.240407458748763e-004 4.240407458748763e-004 wQl4 = wQl14 =

4.240407458748763e-004 4.240407458748763e-004

方法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:4),Y(1:4)), L2=lfun(X(4:7),Y(4:7)),

L3=lfun(X(7:9),Y(7:9)),运行后输出多项式L1，L2，L3.然后输入程

序：

>> syms x

L1=inline('-9070382563990323/56294995342131200*x.^3-7939

56227651389/281474976710656000*x.^2+2253151961697736691/2251799813685248000*x');

L2=inline('-41118170146065463/351843720888320000*x.^3-17

00780706742359/21990232555520000*x.^2+735208215754549839/703687441776640000*x-25688208496779769/2814749767106560000');

L3=inline('-131688789655946847/281474976710656000*x.^2+1

04193434945799191/70368744177664000*x-9309740592764125768759/54295819320043765760000');

[Q141,FCNT141] = quad(L1,0, 0.5498,1.e-14);

[Q142,FCNT142] =quad(L2, 0.5498, 1.0210,1.e-14); [Q143,FCNT143] = quad(L3, 1.0210, pi/2,1.e-14);

Q14= Q141+ Q142+ Q143 ,FCNT14= FCNT141+ FCNT142+ FCNT143 [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);

Q4= Q41+ Q42+ Q43, FCNT4= FCNT41+ FCNT42+ FCNT43 [Ql41,FCNTl41] = quadl(L1,0, 0.5498,1.e-14);

[Ql42,FCNTl42] =quadl(L2, 0.5498, 1.0210,1.e-14); [Ql43,FCNTl43] = quadl(L3, 1.0210, pi/2,1.e-4);

Ql4= Ql41+ Ql42+ Ql43, FCNTl4= FCNTl41+ FCNTl42+ FCNTl43 [Ql141,FCNTl141] = quadl(L1,0, 0.5498,1.e-14);

[Ql142,FCNTl142] =quadl(L2, 0.5498, 1.0210,1.e-14); [Ql143,FCNTl143] = quadl(L3, 1.0210, pi/2,1.e-4);

Ql14= Ql141+ Ql142+ Ql143, FCNTl14= FCNTl141+ FCNTl142+

FCNTl143

wQ14=double(abs(1- Q14)), wQ4=double(abs(1- Q4))

wQl4=double(abs(1- Ql4)), wQl14=double(abs(1- Ql14))

运行后屏幕显示如下

Q14 = FCNT14 = Q4 = FCNT4 = 0.99975250938584 39 0.99975250938584 39 Ql4 = FCNTl4 = Ql14 = FCNTl14 = 0.99975250938584 54 0.99975250938584 54 wQ14 = wQ4 =

2.474906141644029e-004 2.474906141644029e-004 wQl4 = wQl14 =

2.474906141642919e-004 2.474906141642919e-004

方法3 三次样条插值多项式计算数值积分. 输入程序

>> a=0; b= pi/2;

X=[0,0.1571, 0.4712,0.5498,0.6283,0.8639,

1.0210,1.0996,1.5708];

149 ．

第九章 数值积分

fs=[0,0.1565,0.4540,0.5225,0.5878,0.7604,0.8526,0.8910,1.0000];

pp=spline(X,fs);

Qsl=quadl(@ppval,a,b,[],[],pp),wQs1=double(abs(1-Qsl)) Qs8=quad8(@ppval,a,b,[],[],pp),wQs8=double(abs(1-Qs8)) Qs=quad(@ppval,a,b,[],[],pp),wQs=double(abs(1-Qs))

运行后屏幕显示结果如下

Qsl = Qs8 = Qs =

1.00012614786269 1.00012598484705 1.00012589347894 wQs1 = wQs8 = wQs = 1.261478626870538e-004 1.259848470498426e-004 1.258934789407640e-004

9.7 反常积分的计算及其MATLAB程序

9.7.1 无穷积分的符号计算及其MATLAB程序

p??5x例9.7.1 讨论反常积分I??dx的敛散性.

1x4?2解 输入程序

>> syms x

F1=int((5*x)/(x^4+2),x,1,+inf) LimF1= double(F1)

F2=int((5*x.^2)/(x^4+2),x,1,+inf), LimF2= double(F2) F3=int((5*x.^3)/(x^4+2),x,1,+inf), LimF3= double(F3) F8=int((5*x.^8)/(x^4+2),x,1,+inf), LimF8= double(F8)

运行后屏幕显示如下

F1 =

5/8*pi*2^(1/2)-5/4*atan(1/2*2^(1/2))*2^(1/2) LimF1 =

1.6888 F2=

5/4*2^(1/4)*pi-5/8*2^(1/4)*log(1-2^(3/4)+2^(1/2))+5/8

*2^(1/4)*log(1+2^(3/4)+2^(1/2))-5/4*2^(1/4)*atan(2^(1/4)+1)-5/4*2^(1/4)*atan(2^(1/4)-1)

LimF2 = F3 = LimF3 = F8 = LimF8 =

3.9734 inf Inf inf Inf

1???x2pdx,(p??2,0.2,0.5,1,4)的敛散性.

解 当p??2,0.2,0.5,1,4时，输入程序

例9.7.2 讨论反常积分

?1>> syms x

Ff2=int(1/x^(-2*2),x, -inf,-1), F02=int(1/x^(2*0.2),x, -inf,-1)

LimF02= double(F02), F05=int(1/x^(0.5*2),x, -inf,-1) F1=int(1/x^(1*2),x, -inf,-1), F4=int(1/x^(4*2),x, -inf,-1)

运行后屏幕显示

Ff2 = F02 = LimF02 = F05 = Inf -(-1)^(3/5)*inf Inf – Infi -inf

F1 = F4 =

1 1/7

例9.7.3 讨论反常积分

1???x2?2x?bdx（b??3）的敛散性.

??解 输入程序如下

>> syms x

F1=int(1/(x^2+2*x+3),x,-inf, +inf), LimF1= double(F1) F2=int(1/(x^2+2*x-3),x, -inf,0), LimF2= double(F2)

150 ．

第九章 数值积分

F3=int(1/(x^2+2*x-3),x, 0,+inf), LimF3= double(F3) F=F2+F3,LimF=LimF2+ LimF3

运行后屏幕显示如下

F1 = LimF1 =

1/2*pi*2^(1/2) 2.221441469079183e+000 F2 = LimF2 = F3 = LimF3 = NaN NaN NaN NaN F = LimF = NaN NaN

011?2d收敛，而dx不存在，发散.为什么呢？看x?2???x2?2x?3???x?2x?3211了被积函数f(x)?2和g(x)?2的图形就明白了.

x?2x?3x?2x?3即

??输入程序

>> subplot(1,2,1),

x =-10: 0.01:10; y=1./(x.^2+2*x+3);

plot(x,y), grid ,title('y=1/(x^2+2x+3)') subplot(1,2,2),

x =-10: 0.01:10; %z=1./(x.^2+2*x-3);

fplot('1./(x.^2+2*x-3)',[-10 10 -2 2]), grid,title(' y=1/(x^2+2x-3)')

运行后屏幕显示图形（略）.

9.7.2 无穷积分的近似计算及其MATLAB程序

(一) 累积求和无穷积分法及其MATLAB程序

累积求和无穷积分法计算无穷积分的近似值的MATLAB主程序

function [k,suj,wugu,Fjj,WC]=wqjfjx(a,Wucha,m) r=a;k=0; wugu=1; syms t

Y=input( '计算没有结束，请输入被积函数

f(x)=');%exp(-t)./(1+t.^4);

su=int(Y,t,a,a+2) suj= double(su)

while ((wugu>Wucha)&(k

Y=input( '计算没有结束，请输入被积函数f(x)='); F1=int(Y,t,r,r1), intF1= double(F1); suj=suj+intF1,wugu=abs(intF1) Fj=int(Y,t,0,+inf)

Fjj= double(Fj), WC=abs(Fjj- suj) end

disp('k是累加次数,suj是利用累积求和无穷积分法计算f(x)从a到正无

穷的无穷积分的近似值,wugu是f(x)在[r(n)，r(n+1)]上的积分的绝对值,Fjj是f(x)从a到正无穷的无穷积分的值,WC=|Fjj-suj|')

??例9.7.4 计算I?n1?610d，使精确度为，并与精确值比较. xx4?0e(1?x)解 取rn?2，根据（9.72）式计算I???k?0rk??rk?11dx.输入下面的程序： x4e(1?x)>> a=0; Wucha =1.e-6;m=100;

[k,suj,wugu,Fjj,WC]=wqjfjx(a,Wucha,m)

151 ．

第九章 数值积分

运行后屏幕显示

计算没有结束，请输入被积函数f(x)= 输入函数：exp(-t)./(1+t.^4); 运行后屏幕显示

Warning: Explicit integral could not be found.

> In C:\\MATLAB6p5\\toolbox\\symbolic\\@sym\\int.m at line 58 In C:\\MATLAB6p5\\work\\wqjfjx.m at line 5 su =

int(exp(-t)/(1+t^4),t = 0 .. 2) suj =

0.62745967782312 k = 1 r = 2 r1 = 4

计算没有结束，请输入被积函数f(x)= 输入函数：exp(-t)./(1+t.^4); 运行后屏幕显示提示，……，如此继续下去，当wugu<10时，屏幕显示如下内容：

k是累加次数,suj是利用累积求和无穷积分法计算f(x)从a到正无穷的无穷

积分的近似值,wugu是f(x)在[r(n)，r(n+1)]上的积分的绝对值,Fjj是f(x)从a到正无穷的无穷积分的值,WC=|Fjj-suj|

k = 3 suj =

0.63047783491711 wugu =

5.598881094648278e-008 Fjj =

0.63047783491850 WC =

1.386668557756821e-012

(二) 无穷区间的截断法

???6例9.7.5 用数值方法计算I????0e

?3x2dx，使精确度为102?12，并与精确值比较.

解 为估计截断部分

??N?e

?3x2dx的值，利用x?Nx(x?N)时

??N?e

?3x2dx?N?e

?3Nxdx?1?3N2e?10?12. 3N将下面的函数保存为名为nm.m的M文件

function y=nm(x)

y=exp(-3*x.^2)./(3*x); 输入程序

>> N=2:5;y=nm(N)

运行后屏幕显示如下

y = 1.0e-005 *

0.10240353922214 0.00000002088365 0.00000000000000 0.00000000000000

由此可见，如果不计算误差传播，取N?3，可满足I?

?30e

?x2dx?10?12.

152 ．

第九章 数值积分

输入程序

>> x=0:0.1:3; y=exp(-3*x.^2); z=trapz(x,y) syms t,

F1=int(exp(-3*t.^2),t,0,+inf),F3=int(exp(-3*t.^2),t,0,3) intF1= double(F1), Wuchaz=abs(intF1-z)

intF3= double(F3), Wucha3=abs(intF1-intF3) 运行后屏幕显示当N?3时，分别用梯形公式和int函数计算

?30e

?x2dx的近似值z和精

确解 F3（或intF3）及其与I的误差Wuchaz和Wucha3如下

z = F1 = intF1 =

0.51166335397311 1/6*3^(1/2)*pi^(1/2) 0.51166335397324

F3 = intF3 =

1/6*erf(3*3^(1/2))*3^(1/2)*pi^(1/2) 0.51166335397314

Wuchaz = Wucha3 = 1.294520046712933e-013

1.025846074753645e-013

例9.7.6 如果取rk?2k?,k?0,1,2,?，根据（9.72）式，当k取何值时，

??I?sinx?610d的近似值的截断误差为，计算其近似值，并与精确值比较. x2?1?x0解 （1）首先求k的值.输入下面的程序

>> k=ceil(sqrt(10^6/(4*pi)))

运行后屏幕显示

k =

283

?6因此，当k?283时,I的近似值的截断误差为10.

（2）然后计算I的近似值，使截断误差为10，并与精确值比较. 编写并保存名为wqjfjx1.m的计算主程序如下

function [k,suj,wugu]=wqjfjx1(a,m) r=a;k=0; wugu=1; syms t

fun=sin(t)./(1+t.^2); su=int(fun,t,a,2*pi); suj= double(su); while (k

k=k+1,r=2*k*pi; r1=2*(k+1)*pi; syms t

F1=int(fun,t,r,r1); intF1= double(F1); suj=suj+intF1, wugu=abs(intF1), end

输入下面的程序

>> a=0,m=283,[k,suj,wugu]=wqjfjx1(a,m) syms t

fun= sin(t)./(1+t.^2); Fj=int(fun,t,0,+inf), Fjj= double(Fj) WC=abs(Fjj- suj)

运行后屏幕显示计算过程和结果如下：

k = suj = wugu =

1 0.64068026831490 - 0.00000000000000i 0.01602703309687 k = suj = wugu =

2 0.64399926393154 - 0.00000000000000i 0.00331899561664 k = suj = wugu =

3 0.64519499777913 - 0.00000000000000i 0.00119573384758 …………………………………………………………

k = suj = wugu =

283 0.64676080872636 - 0.00000000000000 i 2.223370580497967e-009 Fj =

153 ．

?6

第九章 数值积分

-i*sinint(i)*cosh(1)-cosint(i)*sinh(1)+1/2*i*sinh(1)*pi

Fjj = WC =

0.64676112277913 - 0.00000000000000i 3.140527699008189e-007

(三) 高斯-拉盖尔求积公式及其MATLAB程序

计算n次正交多项式——拉盖尔多项式Un(x)和系数公式An,k(k?1,2,?,n)的MATLAB主程序

function [K,UK,AK]=laguerre(n) syms x

K=0,U0=1,K=1,U1=1-x,u1=1;u2=1-x; duk=diff(u2,x); n=n-1;sun1=1; duk2= duk^2;

A1=(sun1^2)/(x*duk2);A1=simple(A1), sun=1;

for k=1:2:n

u3=(1+2*k-x)*u2-(k^2)*u1;K=k+1, UK =simple(u3), duk1=diff(u3,x); sun2=sun1*K;

duk21= duk1^2;

Ak1=(sun1^2)/(x*duk21); AK =simple(Ak1), u4=(1+2*(k+1)-x)*u3-((k+1)^2)*u2;K=k+2, UK =simple(u4), duk=diff(u4,x);

sun3=sun2*K;sun1=sun3; duk2= duk^2;

Ak2=(sun3^2)/(x*duk2); AK =simple(Ak2), u1=u3;u2=u4;

if k>n,return,end; end

K=k+2; UK; AK;

例9.7.7 根据（9.73）式、（9.74）式和（9.75）式,计算1至7次正交多项

式——拉盖尔多项式、节点的横坐标xj(j?1,2,3,4,5,6,7)和系数An,k(k?1,2,?,7)公式及其在xj(j?1,2,3,4,5,6,7)处的值Ak.

解 （1）建立并保存以laguerre.m文件命名的M文件函数.

（2）计算1至7次正交多项式——拉盖尔多项式和系数An,k(k?1,2,?,7)公式.输入程序

>>n=7,[K,U,Ak]=laguerre(n)

运行后输出结果

K = 0 U0 = 1 K = 1 U1 = 1-x A1 = 1/x K = 2 UK =

2-4*x+x^2 AK =

1/x/(-4+2*x)^2 K = 3 UK =

6-18*x+9*x^2-x^3 AK =

154 ．

第九章 数值积分

4/x/(6-6*x+x^2)^2 K = 4 UK =

-16*x^3+72*x^2-96*x+24+x^4 AK =

9/4/x/(-24+36*x-12*x^2+x^3)^2 K = 5 UK =

-200*x^3+600*x^2-600*x+120+25*x^4-x^5 AK =

576/x/(120*x^2-240*x+120-20*x^3+x^4)^2 K = 6 UK =

720-4320*x+5400*x^2-2400*x^3+450*x^4-36*x^5+x^6 AK =

400/x/(-720+1800*x-1200*x^2+300*x^3-30*x^4+x^5)^2 K = 7 UK =

5040-35280*x+52920*x^2-29400*x^3+7350*x^4-882*x^5+49*x^6-x^7 AK =

518400/x/(-15120*x+5040-4200*x^3+12600*x^2-42*x^5+630*x^4+x^6)^2

（3）求1至7次正交多项式——拉盖尔多项式的零点xk，即求积节点. 输入程序

>>syms x y y1=('1-x=0')

y2=('x^2-4*x+2=0')

y3=('-x^3+9*x^2-18*x+6=0')

y4=('x^4-16*x^3+72*x^2-96*x+24=0')

y5=('-x^5+25*x^4-200*x^3+600*x^2-600*x+120=0')

y6=('720-4320*x+5400*x^2-2400*x^3+450*x^4-36*x^5+x^6=0') y7=('5040-35280*x+52920*x^2-29400*x^3+7350*x^4-882*x^5+4

9*x^6-x^7=0')

y11=solve(y1,x), x1=double(y11) y22=solve(y2,x) , x2=double(y22), y33=solve(y3,x) , x3=double(y33), y44=solve(y4,x) , x4=double(y44), y55=solve(y5,x) , x5=double(y55), y66=solve(y6,x) , x6=double(y66), y77=solve(y7,x) , x7=double(y77),

运行后可以得到（9.74）式中的节点的横坐标xk(k?1,2,3,4,5,6,7)(略).

（4）求在节点的横坐标xj(j?1,2,3,4,5,6,7)处的系数Ak(k?1,2,?,7). 输入程序

>> x1=[1],A1 =1./x1

x2=[3.41421356237309,0.58578643762690],A2 =1./x2./(-4+2*x2).^2 x3=[6.28994508293748 0.41577455678348 2.29428036027904], A3=4./x3./(6-6*x3+x3.^2).^2 x4=[9.39507091230113 4.53662029692113 1.74576110115835

0.32254768961939],

A4 =9/4./x4./(-24+36*x4-12*x4.^2+x4.^3).^2 x5=[0.26356031971814 1.41340305910652 3.59642577104072

7.08581000585884 12.64080084427578],

A5 =576./x5./(120-240*x5+120*x5.^2-20*x5.^3+x5.^4).^2 x6=[0.22284660417926 1.18893210167262 2.99273632605931

5.77514356910451 9.83746741838259 15.98287398060170],

A6 =400./x6./(-720+1800*x6-1200*x6.^2+300*x6.^3-30*x6.^4+x6.^5).^2 x7=[0.19304367656036 1.02666489533919 2.56787674495075

4.90035308452648 8.18215344456286 12.73418029179782 19.39572786226254],

A7

=518400./x7./(5040-15120*x7+12600*x7.^2-4200*x7.^3+630*x7.^4-42*x7.^5+x7.^6).^2

155 ．

第九章 数值积分

运行后输出节点的横坐标xj(j?1,2,3,4,5,6,7)和对应的系数AK （K?1,2,?,7)如下

x1 = 1 A1 = 1 x2 =

3.41421356237309 0.58578643762690 A2 =

0.03661165235168 0.21338834764832 x3 =

6.28994508293748 0.41577455678348 2.29428036027904 A3 =

0.01038925650159 0.71109300992917 0.27851773356924 x4 =

9.39507091230113 4.53662029692113 1.74576110115835 0.32254768961939 A4 =

0.00003370591910 0.00243049428219 0.02233866827736 0.03769713152135 x5 =

Columns 1 through 4

0.26356031971814 1.41340305910652 3.59642577104072 7.08581000585884 Column 5

12.64080084427578 A5 =

Columns 1 through 4

0.52175561058281 0.39866681108317 0.07594244968171 0.00361175867992 Column 5

0.00002336997239 x6 =

Columns 1 through 4

0.22284660417926 1.18893210167262 2.99273632605931 5.77514356910451 Columns 5 through 6

9.83746741838259 15.98287398060170 A6 =

Columns 1 through 4

0.01274901872083 0.01158335641034 0.00314926061317 0.00028886659592 Columns 5 through 6

0.00000725047786 0.00000002495966 x7 =

Columns 1 through 4

0.19304367656036 1.02666489533919 2.56787674495075 4.90035308452648 Columns 5 through 7

8.18215344456286 12.73418029179782 19.39572786226254 A7 =

Columns 1 through 4

0.40931895170127 0.42183127786172 0.14712634865750 0.02063351446872 Columns 5 through 7

0.00107401014328 0.00001586546435 0.00000003170315

用高斯—拉盖尔无穷积分公式计算

???0e

?xf(x)dx的数值积分及其截断误差公式

的MATLAB主程序

function [GL,Y, RGn]=GaussL1(fun,X,A) n=length(X);n2=2*n; Y=feval(fun,X); GL=sum(A.*Y); sun=1; su2n=1; su2n1=1; for k=1:n sun=sun*k; end

for k=1:n2

su2n=su2n*k; end

syms M

RGn= ( sun^2)*M/( su2n);

用高斯—拉盖尔无穷积分公式计算MATLAB主程序

156 ．

???0e

?xf(x)dx的数值积分和误差估计的

第九章 数值积分

function [GL,Y,Rn]=GaussL2 (fun,X,A,fun2n) n=length(X);n2=2*n; Y=feval(fun,X); GL=sum(A.*Y); sun=1; su2n=1; su2n1=1; for k=1:n sun=sun*k; end

for k=1:n2

su2n=su2n*k; end

mfun2n1 =max(fun2n); mfun2n=abs(mfun2n1);

Rn= ( sun^2)*mfun2n /( su2n);

例9.7.8 用高斯—拉盖尔无穷积分公式计算

???01dx，取n=7，再根据x2e(5?x)截断误差公式写出误差公式，并将计算结果与精确值进行比较.

解 （1）建立并保存以fun.m文件命名的M文件函数

function y = fun(x)

y = 1./(5+x.^2);

（2）建立并保存用GaussL2.m文件命名的M文件函数. （3）输入程序

>> syms x

fun2n=diff(exp(-x)./(5+x^2),x,10); fun2ns=simple(fun2n)

（4）运行后，再输入程序

>>X=[0.19304367656036,1.02666489533919,2.56787674495075,

4.90035308452648,

8.18215344456286,12.73418029179782,19.39572786226254];

A=[0.40931895170127,0.42183127786172,0.14712634865750,

0.02063351446872,0.00107401014328,0.00001586546435,0.00000003170315];

x=0:0.001:1000;

fun2ns=exp(-x).*(50889062500*x-29608750000*x.^2+89352450

000*x.^5+69284490000*x.^6+9796290000*x.^7-185554687500*x.^3-70161046875*x.^4-6677396250*x.^8-868165200*x.^10-3734325000*x.^9+14759850*x.^12+1748400*x.^14-79569000*x.^11+7568400*x.^13+36675*x.^16+320*x.^18+x.^20+285840*x.^15+3780*x.^17+20*x.^19+7000234375)./(5+x.^2).^11;

[GL,Y, Rn]=GaussL2(@fun,X,A,fun2ns) syms t

fi=int(exp(-t)/(5+t^2),t,0,+inf); Fs= double (fi), wGL =

double (abs(fi-GL)) 运行后屏幕显示如下

GL =

0.16435187219536 Y =

Columns 1 through 4

0.19852039332099 0.16517893285268 0.08625157652455 0.03446676087003

Columns 5 through 7

0.01389899751582 0.00598231575880 0.00262333436453

Rn =

0.06539672874780 Fs =

0.16435041671009 + 0.00000000000000i wGL =

1.455485269162850e-006

157 ．

第九章 数值积分

(四)高斯—埃尔米特求积公式及其MATLAB程序

高斯—埃尔米特无穷数值积分公式计算

?????e

?x2f(x)dx及其截断误差公式的

MATLAB主程序

function [GH,Y,RHn]=GaussH1(fun,X,A)

n=length(X);n2=2*n; Y=feval(fun,X); GH=sum(A.*Y); sun=1; su2n=1; su2n1=1; for k=1:n sun=sun*k; end

for k=1:n2

su2n=su2n*k; end

syms M

RHn = ( sun*sqrt(pi))*M/( 2^n*su2n);

例9.7.9 用高斯—埃尔米特无穷积分公式计算I?1??2????e

?x2sin3xdx，取

n=5，再根据截断误差公式写出误差公式，并将计算结果与精确值进行比较.

解 （1）建立并保存fun.m文件命名的M文件函数

function y = fun (x)

y = sin(3*x)./sqrt(2*pi);

（2）建立并保存GaussH1.m文件命名的M文件函数. （3）根据表9–8 ，输入程序

>> X=[2.0201828705 -2.0201828705 0.9585724646 -0.9585724646

0];

A=[0.0199532421 0.0199532421 0.3936193232 0.3936193232

0.9453087205];

[GH,Y, RHn]=GaussH1(@fun,X,A) syms x

fi1=int(exp(-x^2).*sin(3*x)./(sqrt(2*pi)),x,0,+inf); fi2=int(exp(-x^2).*sin(3*x)./(sqrt(2*pi)),x, -inf,0); fi=fi1+fi2; Fs= double (fi), wGL = double (abs(fi-GH))

运行后屏幕显示取n=5，用高斯—埃尔米特无穷积分公式数值计算I的结果GH、截断误差公式RHn和精确值Fs及其两者的绝对误差wGL依次如下：

GH = 0 Y =

Columns 1 through 4

-0.08808725583239 0.08808725583239 0.10482362982925 -0.10482362982925

Column 5

0

RHn =

2494507032021713/1361888527316837990400*M Fs = 0 wGL = 0

高斯—埃尔米特无穷数值积分公式计算主程序

function [GH,Y, RHn]=GaussH2(fun,X,A,fun2n) n=length(X);n2=2*n; Y=feval(fun,X); GH=sum(A.*Y); sun=1; su2n=1; su2n1=1; for k=1:n sun=sun*k; end

158 ．

?????e

?x2f(x)dx及其误差估计的MATLAB

第九章 数值积分

for k=1:n2 su2n=su2n*k; end

mfun2n =max(fun2n);

RHn = ( sun*sqrt(pi))* mfun2n /( 2^n*su2n);

例9.7.10 用高斯—埃尔米特无穷积分公式计算

?????x21e(5?x)2dx，取n=5，再

根据截断误差公式写出误差公式，并将计算结果与精确值进行比较.

解 （1）建立并保存fun.m文件命名的M文件函数

function y = fun(x) y = 1./(5+x.^2);

（2）建立并保存GaussH2.m文件命名的M文件函数. （3）输入程序

>> syms x

fun2n=diff(exp(-x.^2)./(5+x.^2),x,10);

fun2ns=simple(fun2n) （4）运行后再输入程序

>> X=[2.0201828705 -2.0201828705 0.9585724646 -0.9585724646

0];

A=[0.0199532421 0.0199532421 0.3936193232 0.3936193232

0.9453087205];

x=-1000:0.001:1000;

fun2ns=32*exp(-x.^2).*(425471484375*x.^2+26348013150*x.^

10-35644522500*x.^8-290495953125*x.^4-262342473750*x.^6+13793263150*x.^12+2639594700*x.^14+34752150*x.^16-88200525*x.^18-18517905*x.^20-1466910*x.^22+32520*x.^24+15840*x.^26+1200*x.^28+32*x.^30-32395781250)./(5+x.^2).^11;

[GH,Y, RHn]=GaussH2(@fun,X,A, fun2ns) syms t

fi1=int(exp(-t.^2)./(5+t.^2),t,0,+inf);

fi2=int(exp(-t.^2)./(5+t.^2),t, -inf,0);

fi=fi1+fi2; Fs= double (fi), wGL = double (abs(fi-GH))

运行后屏幕显示如下

GH =

0.32646126745168 Y =

Columns 1 through 4

0.11011834734512 0.11011834734512 0.16895142009535 0.16895142009535

Column 5

0.20000000000000 RHn =

0.02897468774679 Fs =

0.32640983502878 wGL =

5.143242289669318e-005

9.7.3 无界函数反常积分的符号计算及其MATLAB程序

11例9.7.11 求?dx.

021?x解 （1）输入程序

>> syms x

F=limit(1/(sqrt(1-x^2)),x,1, 'left')

159 ．

第九章 数值积分

运行后屏幕显示结果为

F =

inf 即当x?1?时，被积函数

1.

1?x2??（2）输入程序

>> syms x

F1=int(1/(sqrt(1-x^2)),x,0,1), LimF1= double(F1)

运行后屏幕显示结果为

F1 = LimF1 =

1/2*pi 1.5708

例9.7.12 求?13x5?8?1x2dx.

解 （1）因为被积函数1x2在[?1,1]上除x?0外连续，输入程序

>> syms x

F=limit((3*x^5-8)/ (x^2),x,0)

运行后屏幕显示结果为

F =

-inf

（2）输入程序

>> syms x

F1=int((3*x^5-8)/ (x^2),x,-1,0),

F2=int((3*x^5-8)/ (x^2),x,0,1), F=F1+F2

运行后屏幕显示结果为

F1 = F2 = F =

-inf 35/4 -inf

例9.7.13 讨论反常积分?110xqdx(q?0.5,1,2)的敛散性.

解 （1）因为被积函数

1xq在(0,1]上连续，在x?0无定义，输入程序 >> syms x

LF05=limit(1/(x^0.5),x,0, 'right'), LF1=limit(1/x,x,0, 'right')

LF2=limit(1/(x^2),x,0, 'right')

运行后屏幕显示结果为

LF05 = LF1 = LF2 =

Inf inf inf

（2）输入程序

>> syms x

F05=int(1/(x^0.5),x,0,1),

F1=int(1/x,x,0,1), F2=int(1/(x^2),x,0,1)

运行后屏幕显示结果为

F05 = F1 = F2 =

2 inf inf

9.7.4 无界函数反常积分的近似计算及其MATLAB程序

(一) 无界函数的反常积分的“挖去”法及其MATLAB程序

1 例9.7.14 近似计算I=

?2xlnx01?x2dx，误差在10?6左右.

160 ．

第九章 数值积分

解 输入程序

>> syms r

r=solve('r*(1-(log(r)))=10^(-8)',r), R= double(r)

运行后屏幕显示r的值为

r =

[exp(lambertw(-1/100000000*exp(-1))+1)]

[ exp(lambertw(-1,-1/100000000*exp(-1))+1)] R =

2.71828181845905 0.00000000044374 输入程序

>> r=[2.71828181845905, 0.00000000044374]; y= r.*(1-(log(r)))

运行后屏幕显示r处的函数值为

y =

1.0e-007 *

0.09999995049004 0.10000028042507

取r=0.000 000 000 443 74时，使

12x?6?lnx?10d，所以只需用数值方法x?201?xr计算

?r2xlnx1?x2dx.用梯形公式计算I的近似值Q，输入程序

>> x=0.00000000044374:0.0001:1;

y=(2*x).*abs(log(x))./(1+x.^2);z=trapz(x,y)

运行后屏幕显示

z =

0.41123349222004

用洛巴托求积公式计算，在MATLAB工作区输入程序

>> syms x

Y=inline('(2*x).*abs(log(x))./(1+x.^2)'); Y1=(2*x).*abs(log(x))./(1+x.^2);

[Qql,FCNTl]=quadl(Y, 0.00000000044374,1),

Q1=int(Y1, 0.00000000044374,1), Q1d=double(Q1), juewuQql= Q1d-Qql,

运行后屏幕显示

Qql =

0.41123345753925 FCNTl = 78 Q1d =

0.41123351671206 juewuQql =

5.917280476719355e-008

(二)高斯—切比雪夫求积公式及其MATLAB程序

高斯—切比雪夫求积公式计算I?序

function [GC,X,Y,RCn]=GaussC1(fun,n)

n2=2*n; X=zeros(1,n); A=(pi/n)*ones(1,n); for k=1:n

X(k)=-cos((2*k-1)*pi/(2*n)); end

Y=feval(fun,X); GC=sum(A.*Y); su2n=1; for k=1:n2

su2n=su2n*k;

161 ．

?1?1f(x)1?x2dx及其截断误差公式的MATLAB主程

第九章 数值积分

end

su2n; syms df2n,RCn = (2*pi)*df2n/( 2^(2*n)*su2n); 例9.7.16 用高斯—切比雪夫求积公式计算I??1?12?x1?x2dx，取n=15，再根据

截断误差公式写出误差公式，并将计算结果与精确值进行比较.

解 输入程序

>> n=15; [GC,X,Y,RCn]=GaussC1(@fun,n) syms x

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

运行后屏幕显示如下

GC =

4.36887628549240 X =

Columns 1 through 9

-0.99452189536827 -0.95105651629515 -0.86602540378444 -0.74314482547739 -0.58778525229247 -0.40673664307580 -0.20791169081776 -0.00000000000000 0.20791169081776

Columns 10 through 15

0.40673664307580 0.58778525229247 0.74314482547739 0.86602540378444 0.95105651629515 0.99452189536827

Y =

Columns 1 through 9

1.00273531135177 1.02417941968429 1.06488243304863 1.12109552426303 1.18836641979969 1.26224536320170 1.33868902631726 1.41421356237309 1.48590433434248

Columns 10 through 15

1.55136605708511 1.60865945814907 1.65624419258677 1.69293396320838 1.71786393998336 1.73046869239760

RCn =

1/142406544757979148169424935503213094764544*pi*df2n Fs =

4.36887628549240 wGL =

8.919048090161604e-016

n个节点的高斯—切比雪夫求积公式计算程序

?1?1f(x)1?x2dx及其误差估计的MATLAB主

function [GC, X, Y, RCn]=GaussC2(fun,n ,df2n) n2=2*n; X=zeros(1,n); A=(pi/n)*ones(1,n); for k=1:n

X(k)=-cos((2*k-1)*pi/(2*n)); end

Y=feval(fun,X); GC=sum(A.*Y); su2n=1; for k=1:n2

su2n=su2n*k; end

su2n; mfun2n =max(df2n);

RCn = (2*pi)*mfun2n /( 2^(2*n)*su2n);

例9.7.17 用高斯—切比雪夫求积公式计算I??1?1e?x21?x2dx，取n=6，估计其

误差，并将计算结果与精确值进行比较.

解 输入程序

>> x=-1:0.0001:1;

df2n=64*exp(-x.^2).*(10395-124740*x.^2+207900*x.^4-11088

0*x.^6+23760*x.^8-2112*x.^10+64*x.^12);

n=6; [GC, X, Y, RCn]=GaussC2 (@fun,n ,df2n)

162 ．

第九章 数值积分

syms t a b fi1=int(exp(-t.^2)./sqrt(1-t.^2),t,0,a); fia=limit(fi1,a,1)

fi2=int(exp(-t.^2)./sqrt(1-t.^2),t, -b,0); fib=limit(fi2,b,1), fi=fia+fib;

Fs= double (fi), wGL = double (abs(fi- GC))

运行后屏幕显示结果依次如下

GC =

2.02643676313532 X =

Columns 1 through 4

-0.96592582628907 -0.70710678118655 -0.25881904510252 0.25881904510252 Columns 5 through 6

0.70710678118655 0.96592582628907 Y =

Columns 1 through 4

0.39336682642326 0.60653065971263 0.93520708016087 0.93520708016087 Columns 5 through 6

0.60653065971263 0.39336682642326 RCn =

2.130528872063391e-006

Warning: Explicit integral could not be found. fia =

1/2*pi*hypergeom([1/2],[1],-1)

Warning: Explicit integral could not be found. fib =

int(exp(-t^2)/(1-t^2)^(1/2),t = -1 .. 0) Fs = wGL =

2.02643806694936 1.303814033355041e-006

9.8 多重积分的计算及其MATLAB程序

9.8.1 二重积分的符号计算及其MATLAB程序

例9.8.1 计算??cos(x?y)d?，其中Dxy是由曲线2xy?1,y?2x,x?2.5所围

Dxy成的平面区域.

解 （1）画出积分区域的草图.输入程序

>>x=0.001:0.001:3;

y1=1./(2*x);y2=sqrt(2*x); y=-0.5:0.01:3;

plot(x,y1,'b-',x,y2,'M',2.5,y,'g-'), axis([-0.5 3 -0.5 3]) title('由y1=1/2x,y2=sqrt(2x)和 x=2.5 所围成的积分区域Dxy')

运行后屏幕显示图形（略）.

（2）确定积分限.输入程序

>> syms x y

y1=('2*x*y=1'); y2=('y-sqrt(2*x)=0'); [x,y]=solve(y1,y2,x,y) 运行后屏幕显示两条曲线2xy?1,y?2x的交点如下

x =1/2，y =1 （3）输入计算程序

>> syms x y

f=cos(x+y);y1=1/(2*x);y2=sqrt(2*x); jfy=int(f,y,y1,y2); jfx=int(jfy,x,0.5,2.5); jf2=double(jfx)

运行后屏幕显示如下

Warning: Explicit integral could not be found.

> In D:\\MATLAB6P5\\toolbox\\symbolic\\@sym\\int.m at line 58 jf2 =

-1.83209375329577

163 ．

第九章 数值积分

因此，所求的

Dxy??cos(x?y)d?的近似值为-1.832 093 753 295 77.

sin(x?y)d?，其中Dxy是由曲线x?y2,y?x?2所围成的平??x?yDxy例9.8.2 计算

面区域.

解 （1）画出积分区域的草图.输入程序

>>syms x y

f1=x-y^2;f2=x-y-2; ezplot(f1),hold on

ezplot(f2),hold off,axis([-0.5 5 -1.5 3]) title('由x=y^2和 y=x-2 所围成的积分区域Dxy')

运行后屏幕显示图形（略）.

（2）确定积分限.输入程序

>> syms x y

y1=('x-y^2=0'); y2=('x-y-2=0'); [x,y]=solve(y1,y2,x,y)

运行后屏幕显示两条曲线x?y2,y?x?2的交点如下

x = y =

[ 1] [ -1] [ 4] [ 2]

（3）输入计算程序

>> syms x y

f=sin(x+y)/(x+y);x1=y^2;x2=y+2; jfx=int(f,x,x1,x2); jfy=int(jfx,y,-1,2); jf2=double(jfy)

运行后屏幕显示如下

Warning: Explicit integral could not be found. jf2 =

1.97124962844910

因此，所求的

sin(x?y)d??1.971 249 628 449 10. ??x?yDxy

9.8.2 二重积分的梯形公式及其MATLAB程序

（一） 二重积分的基本梯形公式及其MATLAB程序 例9.8.3 根据基本梯形公式编写MATLAB程序，计算I?sin(x?y)d?近似??x?yDxy值，并将计算结果与精确值比较.其中Dxy是矩形区域0?x?4,?1?y?2.

解 编写并输入下列MATLAB程序

>> a=0;b=4;c=-1;d=2;

fac=sin(a+c)/(a+c);fad=sin(a+d)/(a+d);

fbc=sin(b+c)/(b+c);fbd=sin(b+d)/(b+d); h=(b-a)*(d-c); F=fac+fad+fbc+fbd; T2=(h*F)/4 syms x y

bjh=sin(x+y)/(x+y); jfx=int(bjh,x,a,b);

jfy=int(jfx,y,c,d); I2=double(jfy), juewu=abs(I2-T2)

运行后屏幕显示如下

T2 = I2 = juewu =

3.88977135362262 3.64403676611252 0.24573458751009

(二) 二重积分的复合梯形公式及其MATLAB程序

e?(x?y)例9.8.4 根据（9.93）式编写MATLAB程序计算I???d?近似值，分别

x?yDxy 164 ．

22第九章 数值积分

取(n,m)?(2,2),(12,12),(112,112),(1112,1112)，并将计算结果与精确值比较.观察节点的个数与误差的关系.问取(n,m)??(111112,111112)?时，计算机能计算吗，为什么？其中Dxy是矩形区域0?x?4,1?y?2.

解 （1）编写下列MATLAB主程序并保存名为r2trapf2.m的M文件

function [T2,I2,juewu]=r2trapf2(a,b,n,c,d,m)

x=zeros(1,n);y=zeros(1,m); P=zeros(n,m); fyx=zeros(n,m); hx=(b-a)/(n-1);hy=(d-c)/(m-1);

x(1)=a;x(n)=b;y(1)=c;y(m)=d;

for k=1:n-1 x(k+1)=a+k*hx;

for j=1:m-1 y(j+1)=c+j*hy;

fyx(k+1,j+1)=exp(-(x(k+1)^2+y(j)^2))/(x(k+1)+y(j)); fyx(1,j+1)=exp(-(x(1)^2+y(j+1)^2))/(x(1)+y(j+1)); fyx(k+1,1)=exp(-(x(k+1)^2+y(1)^2))/(x(k+1)+y(1)); fyx(1,1)=exp(-(x(1)^2+y(1)^2))/(x(1)+y(1)); fyx(1,m)=exp(-(x(1)^2+y(m)^2))/(x(1)+y(m)); fyx(n,1)=exp(-(x(n)^2+y(1)^2))/(x(n)+y(1)); fyx(n,m)=exp(-(x(n)^2+y(m)^2))/(x(n)+y(m)); P(k+1,j+1)=4;P(1,j+1)=2;P(n,j+1)=2; P(k+1,1)=2; P(k+1,m)=2;

P(1,1)=1; P(1,m)=1; P(n,1)=1; P(n,m)=1; end end

P;fyx; F=P.*fyx; F1=sum(F); F2=sum(F1); T2=hx*hy*F2/4; syms t w

bjh=exp(-(w^2+t^2))/(w+t); jfx=int(bjh,w,a,b);

jfy=int(jfx,t,c,d); I2=double(jfy); juewu=abs(I2-T2); （2）输入下列MATLAB程序

>>a=0;b=4;n2=2;c=1;d=2;

m2=2;[T2,I2,juewu2]=r2trapf2(a,b,n2,c,d,m2) n12=12;;m12=12;

[T12,I12,juewu12]=r2trapf2(a,b,n2,c,d,m12) n112=112;;m112=112;

[T112,I2,juewu112]=r2trapf2(a,b,n112,c,d,m112) n1112=1112; m1112=1112;

[T1112,I2,juewu1112]=r2trapf2(a,b,n1112,c,d,m1112)

运行后屏幕显示如下

Warning: Explicit integral could not be found. T2 = I2 = juewu2 =

0.37703726923921 0.06889208904324 0.30814518019597 T12 = juewu12 = T112 =

0.21706671440746 0.14817462536422 0.07079360039872 ,juewu112 = T1112 = ,juewu1112 =

0.00190151135548 0.06908616401326 1.940749700208577e-004

由此可见，节点的个数越多，近似值与精确值的绝对误差越小.但是，取

计算机就不能计算了，因为计算机内存不够.见下例： (n,m)??(111112,111112)?时，

>> a=0;b=4;n2=111112;c=1;d=2;m2=111112; [T2,I2,juewu2]=r2trapf2(a,b,n2,c,d,m2) ??? Error using ==> zeros

Product of dimensions is greater than maximum integer.

Error in ==> D:\\MATLAB6p5\\work\\r2trapf2.m On line 2 ==> x=zeros(1,n);y=zeros(1,m); P=zeros(n,m); fyx=zeros(n,m);

165 ．

第九章 数值积分

用复合梯形公式（9.93）数值计算二重积分I2=

Dxy??f(x,y)d?的MATLAB主程序

function [T2,I2,juewu]=r2trapf2(fun,a,b,n,c,d,m)

x=zeros(1,n);y=zeros(1,m); P=zeros(n,m); fyx=zeros(n,m); hx=(b-a)/(n-1);hy=(d-c)/(m-1);

x(1)=a;x(n)=b;y(1)=c;y(m)=d;

for k=1:n-1 x(k+1)=a+k*hx;

for j=1:m

y(j+1)=c+j*hy; fyx(k+1,j+1)=feval(fun,x(k+1),y(j)); fyx(1,j+1)=feval(fun,x(1),y(j+1)); fyx(k+1,1)=feval(fun,x(k+1),y(1)); fyx(1,1)=feval(fun,x(1),y(1)); fyx(1,m)=feval(fun,x(1),y(m)); fyx(n,1)=feval(fun,x(n),y(1)); fyx(n,m)=feval(fun,x(n),y(m));

P(k+1,j+1)=4; P(1,j+1)=2;P(n,j+1)=2; P(k+1,1)=2; P(k+1,m)=2;

P(1,1)=1; P(1,m)=1; P(n,1)=1; P(n,m)=1; end end

P;fyx; F=P.*fyx; F1=sum(F); F2=sum(F1); T2=hx*hy*F2/4; syms t w

bjh=feval(fun,w,t); jfx=int(bjh,w,a,b); jfy=int(jfx,t,c,d);

I2=double(jfy); juewu=abs(I2-T2);

9.8.3 矩形域上的辛普森公式及其MATLAB程序

调用格式一：Q2=dblquad(FUN,a,b,c,d)

调用格式二：Q2=dblquad(FUN,a,b,c,d，tol)

调用格式三：Q2=dblquad(FUN,a,b,c,d，tol，@QUADL) 调用格式四：Q2=dblquad(FUN,a,b,c,d，tol，@MYQUADF)

调用格式五：Q2=dblquad(FUN,a,b,c,d，tol，@QUADL,P1,P2,...) 调用格式六：Q2=dblquad(FUN,a,b,c,d，[ ]，[ ],P1,P2,...)

与 Q2=dblquad(FUN,a,b,c,d，1.e-6,@QUAD,P1,P2,...) 相同.

例9.8.5 分别用MATLAB函数dblquad的调用格式二和三计算

e?(x?y)?4I???d?的值，取误差限为tol=10，并将计算结果与精确值比较.其中Dxy

x?yDxy是矩形区域0?x?4,1?y?2.

解 建立并保存被积函数的M文件

function z=integrnd(x,y)

z=exp(-(x.^2+y.^2))./(x+y); 在MATLAB工作窗口输入程序

>> a=0;b=4;c=1;d=2; Q2=dblquad(@integrnd,a,b,c,d,1.e-4) QL2=dblquad(@integrnd,a,b,c,d,1.e-4,@quadl) syms t w

bjh=exp(-(w^2+t^2))/(w+t); jfx=int(bjh,w,a,b); jfy=int(jfx,t,c,d); I2=double(jfy)

Juewu2=abs(I2-Q2), JuewuL2=abs(I2-QL2)

运行后屏幕显示如下

Q2 = QL2 =

0.0689 0.0689

Warning: Explicit integral could not be found.

166 ．

22第九章 数值积分

I2 = Juewu2 = JuewuL2 = 0.0689 7.8537e-006 4.5212e-006

例9.8.6 分别用MATLAB函数dblquad的调用格式二和三计算

I???sin(x?y)?1?4d?的值，分别取误差限为tol=10，10，并将计算结果与精确值

x?yDxy比较.其中Dxy是矩形区域0?x?4,?1?y?2. 解 在MATLAB工作窗口输入下列MATLAB程序

>>a=0;b=4;c=-1;d=2;

Q12=dblquad(inline('sin(x+y)./(x+y)'),a,b,c,d,1.e-1)

QL12=dblquad(inline('sin(x+y)./(x+y)'),a,b,c,d,1.e-1,@qu

adl)

Q2=dblquad(inline('sin(x+y)./(x+y)'),a,b,c,d,1.e-4)

QL2=dblquad(inline('sin(x+y)./(x+y)'),a,b,c,d,1.e-4,@qua

dl)

syms t w

bjh= sin(w+t)./(w+t); jfx=int(bjh,w,a,b);

jfy=int(jfx,t,c,d); I2=double(jfy)

Juewu12=abs(I2-Q12), JuewuL12=abs(I2-QL12) Juewu2=abs(I2-Q2), JuewuL2=abs(I2-QL2)

运行后屏幕显示如下

Q12 =

3.64404348383673 QL12 =

3.64403680642884 Q2 =

3.64404348383673 QL2 =

3.64403680642884 I2 =

3.64403676611252 Juewu12 =

6.717724212901288e-006 JuewuL12 =

4.031632050427447e-008 Juewu2 =

6.717724212901288e-006 JuewuL2 =

4.031632050427447e-008

?1?4由此可见，分别取误差限为tol=10和10时，计算结果相同，但是用MATLAB函数dblquad的调用格式三比二计算结果误差小.

9.8.4 一般域上二重积分的数值计算及其MATLAB程序

例9.8.7 用MATLAB函数dblquad求直径为8的半球的体积V球，误差为10. 解 在MATLAB工作窗口输入下列MATLAB程序

>> a=-4;b=4;c=-4;d=4;

V1=dblquad(inline('sqrt(max(4^2-(x.^2+y.^2),0))'), a,b,c,d,1.e-4,@quadl)

V=dblquad(inline('sqrt(4^2-(x.^2+y.^2)).*(x.^2+y.^2<=4^2

)'),a,b,c,d,1.e-4)

syms t w

bjh=sqrt(4^2-(w.^2+t.^2));

y1=-sqrt(4^2-t.^2); y2=sqrt(4^2-t.^2);

jfx=int(bjh,w,y1,y2); jfy=int(jfx,t,a,b); I2=double(jfy), JuewuL1=abs(I2-V1)

167 ．

?4第九章 数值积分

Juewu1=abs(I2-V), ezplot('x^2 + y^2 - 16',[-5,5]); axis equal

运行后屏幕显示如下

V1 = V =

1.340434607608882e+002 1.340477821376317e+002 Warning: Explicit integral could not be found. I2 = JuewuL1 =

1.340412865531645e+002 0.00217420772367 Juewu1 =

0.00649558446713

9.8.5 三重积分的计算及其MATLAB程序

(一) 用MATLAB符号计算三重积分 例9.8.8 计算I(f)????(x?e

Vy?sinz)dxdydz，其中积分区域V是由旋转抛

物面z?8?x2?y2，圆柱x2?y2?4和z?0所围成的空间闭区域.

解 （1）画出积分区域V的草图.输入程序

>> [x,y]=meshgrid(-2:0.01:2); z1=8-(x.^2+ y.^2); figure(1)

meshc(x,y,z1) hold on

x=-2:0.01:2; r=2;

[x,y,z]=cylinder(r,30) mesh(x,y,z) hold off

title('由旋转抛物面z=8-(x^2+y^2)，圆柱面 x^2+ y^2=4和z=0 所围

成的积分区域V')

figure(2)

contour(x,y,z,10)

title('由z=8-(x^2+y^2) ，圆柱面 x^2+ y^2=4和z=0 所围成区域V

在x0y面上的投影区域Dxy') 运行后屏幕显示图形（略）.

（2）确定积分限.输入程序

>> syms x y z

f1=('z=8-(x^2+y^2)'); f2=('x^2+ y^2=4');

[x,y,z]=solve(f1,f2,x,y,z)

2222运行后屏幕显示旋转抛物面z?8?x?y和圆柱面x?y?4的交线如下

x = y = z =

[ (4-y^2)^(1/2)] [ y] [ 4] [ -(4-y^2)^(1/2)] [ y] [ 4] （3）输入计算程序

>> syms x y z

f=x+exp(y)+sin(z); z1=0;z2=8-(x^2+y^2); x1=-sqrt(4-y^2);x2=sqrt(4-y^2);

jfz=int(f,z,z1,z2); jfx=int(jfz,x,x1,x2); jfy=int(jfx,y,-2,2);jf2=double(jfy)

运行后屏幕显示如下

Warning: Explicit integral could not be found.

> In C:\\MATLAB6p5p1\\toolbox\\symbolic\\@sym\\int.m at line 58 jf2 =

1.216650998803250e+002

因此，所求的I(f)的近似值为121.665 1.

(二) 用MATLAB数值计算三重积分

168 ．

第九章 数值积分

调用格式一：Q3=triplequad(FUN,a,b,c,d,p,q) 调用格式二：Q3=triplequad(FUN,a,b,c,d,p,q,tol)

调用格式三：Q3=triplequad(FUN,a,b,c,d,p,q,tol,@QUADL) 调用格式四：Q3=triplequad(FUN,a,b,c,d,p,q,tol,@MYQUADF)

调用格式五：Q3=triplequad(FUN,a,b,c,d,p,q,tol,@QUADL,P1,P2,...) 调用格式六：Q3=triplequad(FUN,a,b,c,d,p,q,[ ]，[ ],P1,P2,...)

与 Q3=triplequad(FUN,a,b,c,d,p,q,1.e-6,@QUAD,P1,P2,...) 相同.

例9.8.9 分别用MATLAB函数triplequad 的调用格式二和三计算

I(f)????(x?ey?sinz)dxdydz的值，取误差限为tol=10?4，并将计算结果与精

V确值比较.其中V是三维长方体区域?2?x?2,?2?y?2,0?z?4.

解 建立并保存被积函数的M文件

function u=integrnd1(x,y,z) u=x+exp(y)+sin(z); 在MATLAB工作窗口输入程序

>> a= -2;b=2;c=-2;d=2;p=0;q=4;

Q3=triplequad(@integrnd1,a,b,c,d,p,q,1.e-4)

QL3=triplequad(@integrnd1,a,b,c,d,p,q,1.e-4,@quadl) syms x y z

f=x+exp(y)+sin(z);

jfz=int(f,z,p,q); jfy=int(jfz,y,c,d); jfx=int(jfy,x,a,b); I3=double(jfx), Juewu3=abs(I3-Q3) JuewuL3=abs(I3-QL3)

运行后屏幕显示如下

Q3 =

1.425178451284647e+002,

QL3 =

1.425178451284647e+002 I3 =

1.425178309849224e+002,

Juewu3 =

1.414354233020276e-005 JuewuL3 =

1.414354233020276e-005

169 ．

第九章 数值积分

调用格式一：Q3=triplequad(FUN,a,b,c,d,p,q) 调用格式二：Q3=triplequad(FUN,a,b,c,d,p,q,tol)

调用格式三：Q3=triplequad(FUN,a,b,c,d,p,q,tol,@QUADL) 调用格式四：Q3=triplequad(FUN,a,b,c,d,p,q,tol,@MYQUADF)

调用格式五：Q3=triplequad(FUN,a,b,c,d,p,q,tol,@QUADL,P1,P2,...) 调用格式六：Q3=triplequad(FUN,a,b,c,d,p,q,[ ]，[ ],P1,P2,...)

与 Q3=triplequad(FUN,a,b,c,d,p,q,1.e-6,@QUAD,P1,P2,...) 相同.

例9.8.9 分别用MATLAB函数triplequad 的调用格式二和三计算

I(f)????(x?ey?sinz)dxdydz的值，取误差限为tol=10?4，并将计算结果与精

V确值比较.其中V是三维长方体区域?2?x?2,?2?y?2,0?z?4.

解 建立并保存被积函数的M文件

function u=integrnd1(x,y,z) u=x+exp(y)+sin(z); 在MATLAB工作窗口输入程序

>> a= -2;b=2;c=-2;d=2;p=0;q=4;

Q3=triplequad(@integrnd1,a,b,c,d,p,q,1.e-4)

QL3=triplequad(@integrnd1,a,b,c,d,p,q,1.e-4,@quadl) syms x y z

f=x+exp(y)+sin(z);

jfz=int(f,z,p,q); jfy=int(jfz,y,c,d); jfx=int(jfy,x,a,b); I3=double(jfx), Juewu3=abs(I3-Q3) JuewuL3=abs(I3-QL3)

运行后屏幕显示如下

Q3 =

1.425178451284647e+002,

QL3 =

1.425178451284647e+002 I3 =

1.425178309849224e+002,

Juewu3 =

1.414354233020276e-005 JuewuL3 =

1.414354233020276e-005

169 ．

本文来源：https://www.bwwdw.com/article/ex3o.html