回 帖 发 新 帖 刷新版面

主题:[原创]Matlab求解积分问题探讨

论坛上常看到有人问及积分方面问题,类型相似,特总结如下方法希望对各位学习有所帮助.有错误之处还望各位不吝提出.

一.相关函数:
%符号积分
int(f,v)
int(f,v,a,b)
%数值积分
trapz(x,y) %梯形法沿列方向求函数Y关于自变量X的积分
cumtrapz(x,y) %梯形法沿列方向求函数Y关于自变量X的累计积分
quad(fun,a,b,tol) %采用递推自适应Simpson法计算积分
quad1(fun,a,b,tol) %采用递推自适应Lobatto法求数值积分 
dbquad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol) %二重(闭型)数值积分指令
triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol) %三重(闭型)数值积分指令

二.示例:
例1:计算f(t)=exp(-t^2)在[0,1]上的定积分
本例演示:计算定积分常用方法
>> syms x
int(exp(-x^2),0,1)
 
ans =
 
1/2*erf(1)*pi^(1/2)   %erf为误差函数
>> vpa(int(exp(-x^2),0,1))
 
ans =
 
.7468241328124270
>> d=0.001;x=0:d:1;d*trapz(exp(-x.^2))

ans =

    0.7468

>> quad('exp(-x.^2)',0,1,1e-8)

ans =

    0.7468

例2:计算f(t)=1/log(t)在[0,x],0<x<1上的积分
注意:被积函数于x=0无义,在x-->1^-处为负无穷
本例演示:用特殊函数表示的积分结果,如何用mfun指令
(1)
syms t x
ft=1/log(t);
sx=int(ft,t,0,x)  
sx =
-Ei(1,-log(x))  %完全椭圆函数

(2)
x=0.5:0.1:0.9
sx_n=-mfun('Ei',1,-log(x))              
x =
    0.5000    0.6000    0.7000    0.8000    0.9000
sx_n =
   -0.3787   -0.5469   -0.7809   -1.1340   -1.7758  
(3)%图示被函数和积分函数
clf
ezplot('1/log(t)',[0.1,0.9])            
grid on
hold on
plot(x,sx_n,'LineWidth',3)                
Char1='1/ln(t)';
Char2='{\int_0^x} 1/ln(t) dt';        
title([Char1,'   and    ',Char2])    
legend(Char1,Char2,'Location','SouthWest')  

例3:计算f(t)=exp(-sin(t))在[0,4]上的定积分
注意:本题被函数之原函数无"封闭解析表达式",符号计算无法解题!
本例演示:符号计算有限性
(1)符号计算解法
syms t x
ft=exp(-sin(t))
sx=int(ft,t,0,4)  
ft =exp(-sin(t))

Warning: Explicit integral could not be found.

> In sym.int at 58
sx =
int(exp(-sin(t)),t = 0 .. 4)  

(2)数值计算解法
dt=0.05;                    %采样间隔            
t=0:dt:4;                       %数值计算适合于有限区间上,取有限个采样点               
Ft=exp(-sin(t));        
Sx=dt*cumtrapz(Ft);            %计算区间内曲线下图形面积,为小矩形面积累加得
Sx(end)                %所求定积分值
                                %图示
plot(t,Ft,'*r','MarkerSize',4)
hold on
plot(t,Sx,'.k','MarkerSize',15)
hold off
xlabel('x')
legend('Ft','Sx') 

>>ans =
 3.0632

例4:绘制积分图形,y=2/3*exp(-t/2)*cos(sqrt(3)/2*t);积分s(x)=int(y,t,0,x)于[0,4*pi]上
syms t tao
y=2/3*exp(-t/2)*cos(sqrt(3)/2*t);    
s=subs(int(y,t,0,tao),tao,t);     %获得积分函数            
subplot(2,1,1)                            
                                            %
ezplot(y,[0,4*pi]),ylim([-0.2,0.7])  %单变量符号函数可视化,多变量用ezsurf
grid on                                    
subplot(2,1,2)                            
ezplot(s,[0,4*pi])
grid on
title('s = \inty(t)dt')                    

-------------------- * End * Edited by Guassfans 2007.09.14------------------

回复列表 (共34个回复)

沙发


对函数y=exp(-(x.^2+x+1)/(1+x))从0到正无穷积分,我试了您的几种方法,好像都会报错的,是不是有些函数没法解的?希望能帮忙看看这个能算积分吗?

板凳

支持搂主的倾情奉献,使偶明白了符号积分

3 楼

楼主和各位仁兄,能不能帮我看看我上面的积分为什么不能积啊?有什么方法解决吗?

4 楼

1楼的问题是似无该函数的原函数。
但数值积分肯定是可以的。
0-1: int=.3067601686
0-2: int=.4599633159
0-5: int=.5583068217
0-10: int=.5640928975
0-100: int=.5641346055
0-1000: int=.5641346055

5 楼

不错 顶起

6 楼

四楼能不程序写出来吗?

7 楼

function fy=fxy(x)
fy=exp(-(x.^2+x+1)./(1+x)); % save fxy.m
% main --------
clear,clc
h=.001;p=0;a=0;
while a<100
    p=p+(fxy(a)+fxy(a+h))*h/2;
    a=a+h;
end
p=vpa(p,10)

8 楼

仁兄,写的深奥啊· 小子我实在是看不懂啊#看来我还是多学几年了,再来吧·

9 楼


 perfect!!!

10 楼


无穷求积只能用逼近,主要在于除掉奇点,可以参考自己的精度要求来选取逼近的程度,像无穷远处,或许取为100,或许取为100000,当然实际应用中不可能需要那么大的区间

我来回复

您尚未登录,请登录后再回复。点此登录或注册