数值积分及其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 ．