MATLAB數值積分法

作者：凱魯嘎吉 - 博客園
http://www.cnblogs.com/kailugaji/

一、實驗目的

許多工程技術和數學研究中要用到定積分，如果無法直接算不出精確值（如含在積分方程中的積分）或計算困難但可用近似值近似時，就用數值積分法方法加以解決。常用的算法有：復化梯形、辛甫生（Simpson）、柯特斯（Cotes）求積法; 龍貝格(Romberg)算法；高斯（Gauss）算法。

二、實驗原理

三、實驗程序

下面給出復化Ｓimpson求積法程序（梯形及柯特斯復化求積分程序可比照編制）：

四、實驗內容

選一可精確算值的定積分，用復化的梯形法及復化Ｓimpson求積法作近似計算，並比較結果。

五、解答

1.(程序)

xps.m:

function y=xps(x)
y=x^(3/2);

復化梯形公式：

   trap.m:

function [T,Y,esp]=trap(a,b,n)
h=(b-a)/n;
T=0;
for i=1:(n-1)
    x=a+h*i;
    T=T+xps(x);
end
T=h*(xps(a)+xps(b))/2+h*T;
syms x
Y=vpa(int(xps(x),x,a,b),8);
esp=abs(Y-T);

復化辛甫生（Simpson）公式：

   simpson.m:

function [SI,Y,esp]=simpson(a,b,m)
%a,b為區間左右端點，xps(x)為求積公式，m*2等分區間長度
h=(b-a)/(2*m);
SI0=xps(a)+xps(b);
SI1=0;
SI2=0;
for i=1:((2*m)-1)
    x=a+i*h;
    if mod(i,2)==0
        SI2=SI2+xps(x);
    else 
        SI1=SI1+xps(x);
    end
end
SI=h*(SI0+4*SI1+2*SI2)/3;
syms x
Y=vpa(int(xps(x),x,a,b),8);
esp=abs(Y-SI);

2.(運算結果)

>> [T,Y,esp]=trap(1,2,8)

T =

    1.8636

 
Y =
 
1.8627417
 
 
esp =
 
0.0008089288247354886607354274019599
>> [SI,Y,esp]=simpson(1,2,8)

SI =

    1.8627

 
Y =
 
1.8627417
 
 
esp =
 
0.000000020499792974248975951923057436943

從計算結果看：復化辛普森公式更精確。

3.(拓展（方法改進、體會等）)

MATLAB中有一些內置函數，用於實施自適應求積分，都是根據Gander和Gautschi構造的算法編寫的。

quad：使用辛普森求積，對於低精度或者不光滑函數效率更高

quadl：該函數使用了稱為洛巴托求積（Lobatto Quadrature）的算法，對於高精度和光滑函數效率更高

使用方法：

I=quad(func,a,b,tol);

func是被積函數，a，b是積分限，tol是期望的絕對誤差（如果不提供，默認為1e-6）

例如對於函數f=xe^x在[0,3]上求積分，顯然可以通過解析解知道結果是2e^3+1=41.171073846375336

先創建一個M文件xex.m

內容如下：

function f=xex(x)

f=x.*exp(x);

然后調用：

>> format long

>> format compact

>> quad(@xex,0,3)

ans =

  41.171073850902332

可見有9位有效數字，精度還是蠻高的。

如果使用quadl：

>> quadl(@xex,0,3)

ans =

  41.171074668001779

反而只有7位有效數字