1. 程式人生 > >matlab數值積分中函式積分的4種方法

matlab數值積分中函式積分的4種方法

1.       採用inline行內函數

Matlab中可以有采用幾種不同的方式來指定被積函式。對於簡單的、長度不超過一行的公式採用inline命令比較方便。例如,

可用下面的語句進行計算

>> f=inline('1/sqrt(1+x^4)')    %採用inline行內函數

f =

     行內函數:

     f(x) = 1/sqrt(1+x^4)

>> q=quadtx(f,0,1)

q =

0.9270

從matlab第七版開始,內聯(inline)物件被一種功能更強大的結構匿名函式(anonymousfunction)所替代。在matlab第七版內聯物件還允許被使用,但推薦用匿名函式,因為後者可以生成更高效率的程式程式碼。採用匿名函式,上面的例子變為

>> [email protected](x) 1/sqrt(1+x^4)   %採用匿名函式@(x)

f =

    @(x)1/sqrt(1+x^4)

>> q=quadtx(f,0,1)

q =

0.9270

2.      特殊點不可積函式,採用realmin

如果我們想要計算

可能使用下面的語句

>> f=inline('sin(x)/x')

f =

     行內函數:

     f(x) = sin(x)/x

>> q=quadtx(f,0,pi)

已達到最大遞迴限制 500。使用set(0,'RecursionLimit',N) 可更改此限制。請注意,超出可用堆疊空間可能會使 MATLAB 和/或計算機崩潰。

出錯 inlineeval

此時inline函式積分時好像出現了問題,因為在計算f(0)時,出現了除以0的情況,並且最終產生遞迴限制錯誤。一種補救的方法是,將積分的下限由0變為最小的正浮點數,realmin。

>> q=quadtx(f,realmin,pi)

q =

1.8519

這樣就可以避免遞迴錯誤情況的發生。

3.      使用m檔案寫函式

針對上面出現的問題,另一種方法是使用M檔案,而不是行內函數。建立包含下面程式的檔案sinc.m

function f=sinc(x)

%函式sin(x)/x

if x==0

    f=1;

else

    f=sin(x)/x;

end

然後使用函式控制代碼計算積分

>> Q=quadtx(@sinc,0,pi)

Q =

    1.8519

4.      依賴於引數的積分

一個典型的例子是β函式,它定義為

matlab中已經實現了一個現成的β函式,但我們可以以它為例,說明如何處理積分中的引數。建立一個帶三個引數的行內函數

>> F=inline('t^(z-1)*(1-t)^(w-1)','t','z','w')

或者建立一個M檔案:

function f=betaf(t,z,w)

f= t^(z-1)*(1-t)^(w-1)

並將其命名為betaf.m。

就像任何函式一樣,引數的順序是很重要的。定義被積函式時,必須讓積分變數為其第一個引數。然後給出其他引數的值,作為傳遞給quadtx的附加引數。要計算β(8/3,10/3)

應該先設

>> z=8/3;

>> w=10/3;

>> tol=1e-6;

然後執行命令

>> Q=quadtx(F,0,1,tol,z,w)

Q =

0.0348

或Q=quadtx(@betaf,0,1,tol,z,w)