trapz

1

实验四积分

【实验目的】

1．了解积分的基本概念。

2．掌握积分在计算面积、体积等问题中的应用。

3．了解变步长积分和广义积分

4．学习掌握MATLAB软件有关的命令。

【实验内容】

1．用符号积分命令int计算积分xdxxsin2

2．计算数值积分



1

0

2

2

2

4

1

sin

,dx

x

x

dxx

3．计算数值积分





122

)1(

yx

dxdyyx

【实验准备】

1．积分的有关理论

定积分：积分是微分的无限和，函数)(xf在区间],[ba上的积分定义为











n

i

ii

x

b

a

xfdxxfI

i1

0)max(

)(lim)(

其中.,,2,1),,(,,

1110

nixxxxxbxxxa

iiiiiin





从几何意义上

说，对于],[ba上非负函数)(xf，记分值I是曲线)(xfy与直线bxax,及x轴所围

的曲边梯形的面积。有界连续（或几何处处连续）函数的积分总是存在的。

微积分基本定理（Newton-Leibniz公式）：)(xf在],[ba上连续，且

],[),()('baxxfxF，则有

)()()(aFbFdxxfb

a



这个公式表明导数与积分是一对互逆运算，它也提供了求积分的解析方法：为了求)(xf的

定积分，需要找到一个函数)(xF，使)(xF的导数正好是)(xf，我们称)(xF是)(xf的原

函数或不定积分。不定积分的求法有学多数学技巧，常用的有换元积分和分部积分法。从理

2

论上讲，可积函数的原函数总是存在的，但很多被积函数的原函数不能用初等函数表示，也

就是说这些积分不能用解析方法求解，需用数值积分法解决。

在应用问题中，常常是利用微分进行分析，而问题最终归结为微分的和（即积分）。一

些更复杂的问题是含微分的方程，不能直接积分求解。

多元函数的积分称为多重积分。二重积分的定义为







ij

jiji

yx

G

yxfdxdyyxf

ii

),(lim),(

0)max(22



当),(yxf非负时，积分值表示曲顶柱体的体积。二重积分的计算主要是转换为两次单积分

来解决，无论是解析方法还是数值方法，如何实现这种转换，是解决问题的关键。

平面曲线],[),(),(battyytxx的长度

b

a

dttytxL22)(')('

空间曲线],[),(),(),(battzztyytxx的长度

b

a

dttztytxL222)(')(')('

曲面Gyxyxzz),(),,(的面积















G

dxdy

y

z

x

z

S22)()(1

2．积分的数值方法

梯形法：将],[ba划分为若干小区间,.

10

bxxxa

n

则









n

i

x

x

b

a

i

i

dxxfdxxfI

11

)()(

在每一小区间],[

1ii

xx



上)(xf近似为一直线，用弦代替，有

))()((

2

)(

1

1

1

ii

ii

x

x

xfxf

xx

dxxfi

i













从而













n

i

ii

iixfxf

xx

I

1

1

1))()((

2

称为梯形公式。通常将区间],[ban等分，ihax

n

ab

h

i





,，













1

1

1

))(

2

)()(

(

n

i

in

xf

afbf

hTI

可以证明，当n时由上述公式给出的梯形法是收敛的。

3

变步长积分法：以上介绍的梯形法中，划分是给定的。在实际应用中，等分数n往往是

难以确定的。下面介绍变步长梯形法，其思路是对于给定的精度，从1n开始，等分数逐

次加倍,8,4,2n，直至

nn

TT

2

为此，一般说来，这样的做法会浪费计算量，幸运

的是，

n

T

2

与

n

T有如下的递推关系







n

i

nn

ihaf

h

TT

1

2

)

2

1

(

22

1

其中

n

ab

h



。

重积分：重积分的数值计算可通过若干次单积分的组合实现，如对于二重积分



G

dxdyyxfI),(

先化为二次计分

)(

)(

),(xd

xc

b

a

dyyxfdxI

利用梯形法，先将],[ba区间m等分，.,,1,0,,miihax

m

ab

h

xix





利用梯形积分

公式可得









1

1

)(

)(

.),()()),())()((

2

1

(

m

i

xd

xc

iiix

dyyxfxGxGbGaGhIi

i

再将)](),([

ii

xdxc区间n等分，.,,1,0),(,

)()(

)(njijhay

n

xcxd

ih

yij

ii

y





利用

梯形积分公式可得







1

1

)).,()))(,())(,((

2

1

)(()(

n

j

ijiiiiiyi

yxfxdxfxcxfihxG

3．积分的MATLAB命令

MATLAB中主要用int进行符号积分，用trapz,dblquad,quad,quad8等进行数值积分。

int(s)符号表达式s的不定积分

int(s,x)符号表达式s关于变量x的不定积分

int(s,a,b)符号表达式s的定积分，a,b分别为积分的上、下限

int(s,x,a,b)符号表达式s关于变量x的定积分，a,b分别为积分的上、下限

trapz(x,y)梯形积分法，x时表示积分区间的离散化向量，y是与x同维数的向量，

表示被积函数，z返回积分值。

quad8(‘fun’,a,b,tol)变步长数值积分，fun表示被积函数的M函数名，a,b分别为积

分上、下限，tol为精度，缺省至为1e-3.

fblquad(‘fun’,a,b,c,d)矩形区域二重数值积分，fun表示被积函数的M函数名，a,b

分别为x的上、下限，c,d分别为y的上、下限.

可以用helpint,helptrapz,helpquad等查阅有关这些命令的详细信息

4

【实验方法与步骤】

练习1用符号积分命令int计算积分xdxxsin2，MATLAB代码为：

>>clear;symsx;

>>int(x^2*sin(x))

结果为

ans=-x^2*cos(x)+2*cos(x)+2*x*sin(x)

如果用微分命令diff验证积分正确性，MATLAB代码为：

>>clear;symsx;

>>diff(-x^2*cos(x)+2*cos(x)+2*x*sin(x))

结果为

ans=x^2*sin(x)

练习2计算数值积分



1

0

2

2

2

4

1

sin

,dx

x

x

dxx.先用梯形积分法命令trapz计算积分



2

2

4dxx，MATLAB代码为：

>>clear;x=-2:0.1:2;y=x.^4;%积分步长为0.1

>>trapz(x,y)

结果为

ans=12.8533

实际上，积分

2

2

4dxx的精确值为8.12

5

64

。如果取积分步长为0.01,MATLAB代码为：

>>clear;x=-2:0.01:2;y=x.^4;%积分步长为0.01

>>trapz(x,y)

结果为

ans=12.8005

可用不同的步长进行计算，考虑步长和精度之间的关系。一般说来,trapz是最基本的数

值积分方法,精度低,适用于数值函数和光滑性不好的函数.

如果用quad8命令计算积分

2

2

4dxx，得先编写被积函数2)(xxf的M函数

%M函数fun1.m

functiony=fun1(x)

y=x.^4;

保存后,在命令窗口用MATLAB代码:

>>clear;

>>quad8('fun1',-2,2)

>>vpa(quad8('fun1',-2,2),10)%以10位有效数字显示结果

结果为

ans=12.8000

ans=12.80000000

对于变步长数值积分,常用的有quad,quad8两种命令,quad使用自适应步长Simpson法,

5

quad8使用自适应步长8阶Newton-Cotes法,我们建议用quad8,它不但精度较高,且对假收敛

(见习题6)和假奇异积分具有一定的适应性,而quad较差..

如果用符号积分法命令int计算积分

2

2

4dxx，MATLAB代码为：

>>clear;symsx;

>>int(x^4,x,-2,2)

结果为

ans=64/5

对于定积分,



1

0

2

1

sin

dx

x

x

,用步长为0.1,0.01的梯形积分命令trapz计算的结果分别为

(MATLAB代码略去)

ans=0.1811

ans=0.1808

用quad8命令计算结果为(显示10位有效数字)

ans=.1807896039

如果用符号积分法命令int计算积分



1

0

2

1

sin

dx

x

x

，MATLAB代码为：

>>clear;symsx;

>>int(sin(x^2)./(x+1),x,0,1)

>>vpa(int(sin(x^2)./(x+1),x,0,1),10)

结果为:

?Warning:Explicitintegralcouldnotbefound.

>InD:Ronghwmatlabtoolboxsymbolic@e58

InD:e3

ans=int(sin(x^2)/(x+1),x=0..1)%说明求不出解析解

ans=.1807896039%可求近似数值解

这说明int求不出符号解,但可求数值解.一般说来,int的功能虽然很强大,但计算速度慢,

数值计算效率不好.

练习3计算数值积分





122

)1(

yx

dxdyyx,可将此二重积分转化为累次积分









1

1

1

1

1

2

2

22

)1()1(x

x

yx

dyyxdxdyyx

先编写四个M函数文件,

%二重积分算法文件dblquad2.m

functionS=dblquad2(f_name,a,b,c_lo,d_hi,m,n)

6

%其中f_name为被积函数字符串,'c_lo'和'd_hi'是y的下限和上限函数,都是x的标量函数;a,b

分别为x的下限和上限;m,n分别为x和y方向上的等分数(缺省值为100).

ifnargin<7,n=100;end

ifnargin<6,m=100;end

ifm<2|n<2

error('Numnerofintervalsinvalid');

end

mpt=m+1;hx=(b-a)/m;x=a+(0:m)*hx;

fori=1:mpt

ylo=feval(c_lo,x(i));yhi=feval(d_hi,x(i));

hy=(yhi-ylo)/n;

fork=1:n+1y(i,k)=ylo+(k-1)*hy;f(i,k)=feval(f_name,x(i),y(i,k));end

G(i)=trapz(y(i,:),f(i,:));

end

S=trapz(x,G);

%被积函数eg3_fun.m

functionz=eg3_fun(x,y)

z=1+x+y;

%积分下限函数eg3_low.m

functiony=eg3_low(x)

y=-sqrt(1-x^2);

%积分上限函数eg3_up.m

functiony=eg3_up(x)

y=sqrt(1-x^2);

保存后,在命令窗口用MATLAB代码:

>>clear;

>>dblquad2('eg3_fun',-1,1,'eg3_low','eg3_up')

结果为

ans=3.1383

实际上，我们知道此二重积分的精确值为

14159265.3)1(

122







yx

dxdyyx。

为了得到更精确的数值解，需将区间更细化，比如x和y方向等分为1000分，MATLAB代码：

>>clear;dblquad2('eg3_fun',-1,1,'eg3_low','eg3_up',1000,1000)

结果为ans=3.1415。

此题也可用int符号计算求解，MATLAB代码为：

>>clear;symsxy;

>>iy=int(1+x+y,y,-sqrt(1-x^2),sqrt(1-x^2));

>>int(iy,x,-1,1)

结果为

ans=pi

练习4（假奇异积分）计算数值积分

1

1

3

1

dxx，如果用梯形法计算，MATLAB代码为：

>>clear;x=-1:0.01:1;y=x.^(1/3);

7

>>trapz(x,y)

结果为

ans=1.1241+0.6490i

是一个复数。实际上，我们知道0

1

1

3

1



dxx。这个错误是由于数值方法的特点造成的，数

值方法对3

1

x是通过)3/)exp(ln(x计算，当0x时就会出现复数，这种情况称为假奇异积

分。这类情况在适当定义被积函数（当0x时，用)))((3

1

x后可用quad,quad8求解。先

编写M函数

functiony=fun5(x)

y=x.^(1/3);

ifx<0,y=-(-x).^(1/3);end

保存后,在命令窗口用MATLAB代码:

>>clear;quad8('fun5',-1,1)

结果为

ans=5.9940e-004+3.4606e-004i

练习6（广义积分）计算广义积分



dx

x

xI)

50

exp(sin

2

。广义积分数值求解是个

较困难的问题，通常数值积分方法都不适用。若先编一个M函数

functiony=fun6(x)

y=exp(sin(x)-x.^2/50);

再用命令

>>quad8('fun6',-inf,inf)

结果为

ans=NaN

是无穷大，错误。为能正确求解，考虑到











N

N

N

dx

x

xdx

x

xI)

50

exp(sinlim)

50

exp(sin

22

但问题是N去多大？比如取1010N，用命令

>>quad8('fun6',-1e+10,1e+10)

结果为

ans=6.8135e+005

是一个明显错误的大数。正确的方法是，先取一个适当大的N（比如10），计算一个I值；

然后将N以一适当的倍数（如2）增加计算出新的I值，直至前后两次差异小于给定精度为

至。MATLAB代码为：

>>clear;n=10;r=2;eps=1e-5;t0=inf;t1=quad8('fun6',-n,n);

>>whileabs(t0-t1)>eps,t0=t1;n=n*r;t1=quad8('fun6',-n,n),end

结果为15.8678。

也可用符号积分命令int进行计算，MATLAB代码为：

>>symsx;

>>y=int(exp(sin(x)-x^2/50),-inf,inf);

>>vpa(y,10)

结果为15.86778263。

8

练习7（暇积分）计算暇积分



1

0)!(sinxx

dx

，若先编一个M函数

functiony=fun7(x)

y=1./((1+sin(x)).*sqrt(x));

再用命令

>>clear;quad8('fun7',0,1)

结果为

ans=NaN

是无穷大，错误，因为此积分显然收敛。如果用

>>quad8('fun7',1e-5,1)

>>quad8('fun7',1e-10,1)

结果为1.5811和3.2809，系统警告结论不可靠，那个正确？

可用符号积分命令int进行计算，MATLAB代码为：

>>symsx;

>>y=int(1/((1+sin(x))*sqrt(x)),0,1);

>>vpa(y,10)

结果为1.586625183。

【练习与思考】

1．(不定积分)用int计算下列不定积分，并用diff验证

dxxx2sin，

x

dx

cos1

，

1xe

dx

，xdxarcsin，xdx3c

2．(定积分)用trapz,quad8,int计算下列定积分

1

0

sin

dx

x

x

,1

0

dxxx,2

0

)2sin(dxxex,

1

0

2dxex

3．(椭圆的周长)用定积分的方法计算椭圆1

49

22



yx

的周长

4．（二重积分）计算数值积分





yyx

dxdyyx

222

)1(

5．（假奇异积分）用trapz,quad8计算积分

1

1

30cosxdxx。，会出现什么问题？分析原因，

并求出正确的解。

6．（假收敛现象）考虑积分

kdxxkI

0

sin)(，

（1）用解析法求)(kI;

（2）分别用trapz,quad和quad8求)6(),4(II和)8(I,发现什么问题？

7．(Simpson积分法)编制一个定步长Simpson法数值积分程序.计算公式为

)42424(

3114321



nnnn

fffffff

h

SI

其中n为偶数,.1,,2,1),)1((,



nihiaff

n

ab

n

i



8．(广义积分)计算广义积分







dx

x

x

4

2

1

)exp(

,1

0

)tan(

dx

x

x

,



1

0

21

sin

dx

x

x

并验证公式

9









,1

2

)

2

exp(

2

dx

x





02

sin

dx

x

x

.