Matlab 數值積分

Matlab 數值積分
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<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')      </X<1上的積分