MATLAB學習筆記：數值積分

當

（1）被積函式的原函式不能用初等函式表示。

（2）被積函式難以用公式表示，而是用圖形或表格給出的。

就應該建立定積分的近似計算方法：數值積分方法。

梯形法：

z=trapz(x,y)

>> x=0:0.5:1;
>> y=exp(-x.^2);
>> z=trapz(x,y)

z =

    0.7314

>> x=0:0.05:1;
>> y=exp(-x.^2);
>> z=trapz(x,y)

z =

    0.7467

拋物線法：Simpson（辛普森）公式

自適應辛普森公式:

z=quad(f,a,b,tol)

a是積分下限，b是積分上限，tol是計算精度，預設為0.001。

呼叫quad函式時，先要建立一個描述被積函式的函式檔案或語句函式。

以下分別為：行內函數和函式控制程式碼

>> g=inline('exp(-x.^2)');
>> z=quad(g,-1,1)

z =

    1.4936

>> f=@(x) exp(-x.^2);
>> z=quad(f,-1,1)

z =

    1.4936

分別使用左端點的矩形法，右端點的矩形法和梯形法

>> n=100;
>> x=linspace(0,1,n);
>> y=4./(1+x.^2);

>> jxl=sum(y(1:(n-1)))./n

jxl =

    3.1202

>> jxr=sum(y(2:n))./n

jxr =

    3.1002

>> tx=trapz(x,y)

tx =

    3.1416

還有兩個常用的求數值積分的函式：

z=quadl(f,a,b,tol)  自適應複合Lobatto數值積分法

z=quadgk(f,a,c,tol)  自適應複合Gauss-Kronrod數值積分法，適用於高精度和震盪數值積分，以及廣義數值積分

>> g=inline('4./(1+x.^2)');
>> quad(g,0,1)

ans =

    3.1416

>> quadl(g,0,1)

ans =

    3.1416

>> g=@(x) 4./(1+x.^2);
>> quadgk(g,0,1)

ans =

    3.1416

>> f1=@(x) 1./(x.^2);
>> z1=quadgk(f1,1,inf)

z1 =

     1

>> f2=@(x)exp(-x.^2);
>> z2=quadgk(f2,-inf,inf)

z2 =

    1.7725
>> f3=@(x)1./(x.*sqrt(1-log(x).^2));
>> z3=quadgk(f3,1,exp(1))

z3 =

    1.5708