MATLAB 数值积分

作者：未知信息来源:未知 2006-1-26

字体大小:小 中 大网友评论条 进入论坛

8 数值积分5。1 一元函数的数值积分5。1 闭型数值积分【 * 例 5。IS=int(‘exp(-x*x)‘,‘x‘,0,1) % 求解析积分vpa(IS) % 求所得解析积分的 32 位精度近似值IS =1/2*erf(1)*pi^(1/2)ans =。...

5.8 数值积分

5.8.1 一元函数的数值积分

5.8.1.1 闭型数值积分

【 * 例 5.8.1 .1-1 】求，其精确值为。

（1）符号解析法
syms x;IS=int('exp(-x*x)','x',0,1) % 求解析积分
vpa(IS) % 求所得解析积分的 32 位精度近似值
IS =
1/2*erf(1)*pi^(1/2)
ans =
.74682413281242702539946743613185

（2） MATLAB 指令 quad 和 quad8 求积
fun=inline('exp(-x.*x)','x'); % 注意：数组乘符号 .* 的采用必须的。
Isim=quad(fun,0,1),I8=quad8(fun,0,1)
Isim =
0.7468
I8 =
0.7468

（3） 10 数 Gauss 法
Ig=gauss10(fun,0,1)
Ig =
0.7463

（4）样条函数积分法
xx=0:0.1:1.5;ff=exp(-xx.^2); % 产生被积函数的“表格”数据
pp=spline(xx,ff); % 由“表格”数据构成样条函数
int_pp=fnint(pp); % 求样条积分
Ssp=ppval(int_pp,[0,1])*[-1;1] % 据样条函数计算 [0,1] 区间的定积分
Ssp =
0.7468

（5） SIMULINK 积分法

图 5.8.1 .1-1 积分模型 exm5811_1_5.mdl

5.8.1.2 开型数值积分

（1）任意区间上的 10 参数 Gaussian 求积法
[gauss10.m]
function g = gauss10(fun,a,b)
%GAUSS10(fun,a,b) 利用 10 参数 Gauss 求积法近似计算 a<x<b 区间上的积分
% fun 被积函数。字符串、内联函数名、 M 函数文件名都可。
%======================================================
x = [0.1488743390;0.4333953941;0.6974095683; ...
0.8650633667;0.9739065285];
w = [0.2955242247;0.2692667193;0.2190863625; ...
0.1494513492;0.0666713443];
t = .5*(b+a)+.5*(b-a)*[-flipud(x);x];
W = [flipud(w);w];
g = sum(W.*feval(fun,t))*(b-a)/2;

【 * 例 5.8.1 .2-1 】当时，比较解析积分和近似积分。

（1）用符号法求解：
syms x;F=int('cos(x)','x',-1,1),vpa(F)
F =
2*sin(1)
ans =
1.6829419696157930133050046432606

（2）用式 ( 5.8.1 .2-4) 近似计算：
aF=cos(1/sqrt(3))+cos(-1/sqrt(3))
aF =
1.6758

【 * 例 5.8.1 .2-2 】求，准确结果是。

（1）符号法
syms x;IS=int('sqrt(log(1/x))','x',0,1)
Warning: Explicit integral could not be found.
> In D:\MAT53\toolbox\symbolic\@sym\int.m at line 58
In D:\MAT53\toolbox\symbolic\@char\int.m at line 9
IS =
int(log(1/x)^(1/2),x = 0 .. 1)

（2）用 quad 指令求积
ff=inline('sqrt(log(1./x))','x');Isim=quad(ff,0,1)
Warning: Divide by zero.
> In D:\MAT53\toolbox\matlab\funfun\@inline\feval.m at line 30
In D:\MAT53\toolbox\matlab\funfun\quad.m at line 40
Isim =
Inf

（3） 10 参数 Gauss 法
Ig=gauss10(ff,0,1)
Ig =
0.8861

（4）通过变量置换后再进行
令，于是积分成为。
% 对变量置换后的积分利用符号法求解
syms xi;Isxi=int('sqrt(xi)*exp(-xi)','xi',0,inf),vpa(Isxi)
Isxi =
1/2*pi^(1/2)
ans =
.88622692545275801364908374167055

% 对变量置换后的积分利用 quad 求解
gg=inline('sqrt(xi).*exp(-xi)','xi'),Igg=quad(gg,0,inf)
gg =
Inline function:
gg(xi) = sqrt(xi).*exp(-xi)
Igg =
NaN

5.8.2 多重数值积分

5.8.2.2 内积分限为函数的二重积分

[double_int.m]
function SS=double_int(fun,innlow,innhi,outlow,outhi)
%double_int 采用 8 阶多项式近似的二重数值积分
% 内积分变量为 x ，其积分区间上下限可以是变量 y 的表达式。
%fun 被积函数，以 x,y 为积分变量。
%innlow,innhi 内积分区间上下限，可以是标量，或函数文件名字符串
%outlow,outhi 外积分区间上下限，是标量。
y1=outlow;y2=outhi;x1=innlow;x2=innhi;f_p=fun;
SS=quad8( 'G_yi' ,y1,y2,[],[],x1,x2,f_p);

[G_yi.m]
function f=G_yi(y,x1,x2,f_p)
%G_yi 实现 " 内重 " 积分的函数文件，专供 double_int.m 调用。
%y " 外重 " 积分的积分变量
%x1,x2 内积分区间的下、上限。
% 它们由 double_int.m 中 quad8 指令以参数方式传递过来。
%f_p 被积原函数，它由 double_int.m 中 quad8 指令以参数方式传递过来。
y=y(:);n=length(y);
if isstr(x1)==1;xx1=feval(x1,y); else xx1=x1*ones(size(y)); end
if isstr(x2)==1;xx2=feval(x2,y); else xx2=x2*ones(size(y)); end
for i=1:n
f(i)=quad8(f_p,xx1(i),xx2(i),[],[],y(i));
end
f=f(:);

【 * 例 5.8.2 .2-1 】计算。

（1）编写内积分区间上下限的 M 函数文件
[x_low.m]
function f=x_low(y)
f=sqrt(y);

（2）编写被积函数
被积函数函数采用内联函数 ff=inline('x.^2+y.^2','x','y') 表示。

（3）被积函数用内联函数表达时，运行以下指令，即得结果
ff=inline('x.^2+y.^2','x','y'); % 被积函数的内联函数表达
SS=double_int(ff,'x_low',2,1,4)
SS =
9.5810

（4）本题用符号计算很容易获得高精度解
Ssym=vpa(int(int('x^2+y^2','x','sqrt(y)',2),'y',1,4))
Ssym =
9.580952380952381

5.8.3 卷积

5.8.3.5 卷积的MATLAB实现

【 * 例 5.8.3 .5-1 】有序列和。
（A）求两个完整序列的卷积，并图示。（ B ）假设及其后的 4 个非零值未知，而成为截尾序列，求卷积并图示。

% 完整序列卷积
a=ones(1,10);n1=3;n2=12; % 完整 A(n) 序列的非平凡值和区间端点
b=ones(1,8);n3=2;n4=9; % 完整 B(n) 序列的非平凡值和区间端点
c=conv(a,b);nc1=n1+n3;nc2=n2+n4; % 计算卷积和确定卷积非平凡区间端点
kc=nc1:nc2; % 构成非平凡区间的序号自变量

% 截尾序列卷积
aa=a(1:6);nn1=3;nn2=8; % 截尾 A(n) 序列的非平凡值和区间端点
cc=conv(aa,b);ncc1=nn1+n3;
nx=nn2+n4; % “非平凡”区间右端点
ncc2=min(nn1+n4,nn2+n3); % 截尾序列卷积被正确计算区间的右端点
kx=(ncc2+1):nx;kcc=ncc1:ncc2;N=length(kcc);
stem(kcc,cc(1:N),'r','filled')
axis([nc1-2,nc2+2,0,10]),grid,hold on
stem(kc,c,'b'),stem(kx,cc(N+1:end),'g'),hold off

图 5.8.3 .5-1 “完整”序列卷积和“截尾”序列卷积

【 * 例 5.8.3 .5-2 】求函数和的卷积。本例展示：（ A ）符号 Laplace 变换求卷积的理论表示；（ B ） SIMULINK 卷积法的执行过程和它的快速精确性。（ C ）从理论符号解产生相应的理论数值序列。

（1）符号卷积法（得解析结果）
syms tao;t=sym('t','positive'); % 把 t 定义为“取正”符号变量
US1=laplace(exp(-t)); %u(t) 的 L 氏变换
HS1=laplace(t*exp(-t/2)) %h(t) 的 L 氏变换
yt1=simple(ilaplace(US1*HS1)) %L 氏反变换得卷积的理论解
HS1 =

1/(1/2+s)^2
yt1 =
4*exp(-t)+(2*t-4)*exp(-1/2*t)

（2） SIMULINK 卷积法
根据的传递函数，构作 SIMULINK 模型 exm5835_2_2.mdl ，并运行，可从示波器上看到卷积结果。

图 5.8.3 . 5-2-1 卷积解算模型 exm5835_2_2.mdl

（3）比较两种卷积方法的结果，并用图形表示
t=yt2(:,1); %exm5835_2.mdl 示波器所保存的仿真采样时间序列
yyt1=eval(vectorize(char(yt1))); % 理论解的数值序列
[dy,kd]=max(abs(yyt1-yt2(:,2)));
dy12=dy/abs(yyt1(kd))
dy12 =
2.8978e-006

【 * 例 5.8.3 .5-3 】用“零阶”近似法求和的卷积。本例演示：（ A ）连续函数的有限度采样。（ B ）卷积数值计算三个误差（“截尾”误差、“零阶”近似误差、计算机字长误差）的影响。（ C ）卷积“无截尾误差”区间、“非平凡”区间端点的确定。（ D ）离散卷积和连续卷积之间的关系。（ E ）指令 conv 的使用。（ F ）绘图分格线的运用。

（1）离散卷积运算前的准备
由于数值计算只能对有限长度序列进行，所以必须对被卷函数做截尾处理。假如把 5% 作为函数的截断阈值，那么在 t=3 处截尾，在 t=11 处截尾，即取，。

取采样周期。

（2）实施计算
% “零阶”法计算卷积
t2=3;t4=11;T=0.01;
tu=0:T:t2;N2=length(tu);
th=0:T:t4;N4=length(th);
u=exp(-tu);h=th.*exp(-th/2);
tx=0:T:(t2+t4);Nx=length(tx);
yt3=T*conv(u,h);

% 把计算结果与理论值比较
t=tx;yyt1=eval(vectorize(char(yt1)));% 例 5.8.3 .5-2 第（ 1 ）步的计算结果
[dy,kd]=max(abs(yyt1(1:N2)-yt3(1:N2)));
dy13(1)=dy/abs(yyt1(kd));
[dy,kd]=max(abs(yyt1(N2+1:N4)-yt3(N2+1:N4)));
dy13(2)=dy/abs(yyt1(N2+kd));
[dy,kd]=max(abs(yyt1(N4+1:Nx)-yt3(N4+1:Nx)));
dy13(3)=dy/abs(yyt1(N4+kd));

（3）不同区间段的误差及图示
disp(' 与理论结果的相对误差 ')
disp([blanks(4),'[0,3] 段 [3,11] 段 [11,14] 段 ']),disp(dy13)
plot(t,yyt1,':b',t,yt3,'r')
set(gca,'Xtick',[0,3,11,14]),grid

与理论结果的相对误差
[0,3] 段 [3,11] 段 [11,14] 段
0.0068 0.0810 0.6974