一般區域二重、三重積分MATLAB計算方法
阿新 • • 發佈：2019-01-04
這裡討論的計算方法指的是利用現有的MATLAB函式來求解，而不是根據具體的數值計算方法來編寫相應程式。目前最新版的2009a有關於一般區域二重積分的計算函式quad2d，但沒有一般區域三重積分的計算函式，而NIT工具箱似乎也沒有一般區域三重積分的計算函式。
本貼的目的是介紹一種在7.X版本MATLAB（不一定是2009a）裡求解一般區域二重三重積分的思路方法。需要說明的是，在MATLAB的dblquad幫助文件裡已經討論了一種求解一般區域二重三重積分的思路方法，就是將被積函式“延拓”到矩形或者長方體區域，但是這種方法不可避免引入很多乘0運算浪費時間。因此，新的思路將避免這些。由於是呼叫已有的MATLAB函式求解，在求一般區域二重積分時，效率和2009a的quad2d相比有一些差距，但是相對於"延拓"函式的做法，效率大大提高了。下面結合一些簡單例子說明下計算方法。
譬如二元函式f(x,y) = x*y，y從sin(x)積分到cos(x)，x從1積分到2，這個積分可以很容易用符號積分算出結果
如果用的不是2009a，那麼你可以利用NIT工具箱裡的quad2dggen函式。
那麼我們如果既沒有NIT工具箱用的也不是2009a，怎麼辦呢？
答案是我們可以利用兩次quadl函式，注意到quadl函式要求積分表示式必須寫成向量化形式，所以我們構造的函式必須能接受向量輸入。見如下程式碼
x的函式，要能接受向量形式的x輸入，所以構造這個函式的時候考慮到x是向量的情況。
利用arrayfun函式（7.X後的版本都有這個函式，不瞭解這個函式的朋友可以檢視幫助文件，或者百度搜索arrayfun）可以將IntDemo函式精簡成一句程式碼：
首先
而
而
有了這個模板，我們可以方便求其他一般積分割槽域（上下限是函式）形式的二重積分，例如被積函式
f
= @(x,y)
exp(sin(x))*ln(y)，y從5*x積分到x^2，x從10積分到20。
用quad2d,上面介紹的方法，還有dblquad幫助文件裡給的延拓函式的方法
本貼的目的是介紹一種在7.X版本MATLAB（不一定是2009a）裡求解一般區域二重三重積分的思路方法。需要說明的是，在MATLAB的dblquad幫助文件裡已經討論了一種求解一般區域二重三重積分的思路方法，就是將被積函式“延拓”到矩形或者長方體區域，但是這種方法不可避免引入很多乘0運算浪費時間。因此，新的思路將避免這些。由於是呼叫已有的MATLAB函式求解，在求一般區域二重積分時，效率和2009a的quad2d相比有一些差距，但是相對於"延拓"函式的做法，效率大大提高了。下面結合一些簡單例子說明下計算方法。
譬如二元函式f(x,y) = x*y，y從sin(x)積分到cos(x)，x從1積分到2，這個積分可以很容易用符號積分算出結果
-
syms x
-
y
-
int(int(x*y,y,sin(x),cos(x)),1,2) ]
-
結果是
-
-1/2*cos(1)*sin(1)-1/4*cos(1)^2+cos(2)*sin(2)+1/4*cos(2)^2 =
- -0.635412702399943
-
quad2d(@(x,y)
- x.*y,1,2,@(x)sin(x),@(x)cos(x),'AbsTol',1e-12)
如果用的不是2009a，那麼你可以利用NIT工具箱裡的quad2dggen函式。
那麼我們如果既沒有NIT工具箱用的也不是2009a，怎麼辦呢？
答案是我們可以利用兩次quadl函式，注意到quadl函式要求積分表示式必須寫成向量化形式，所以我們構造的函式必須能接受向量輸入。見如下程式碼
-
function
-
IntDemo
-
function f1 = myfun1(x)
-
f1 = zeros(size(x));
-
for k =
-
1:length(x)
-
f1(k) = quadl(@(y)
-
x(k)*y,sin(x(k)),cos(x(k)));
-
end
-
end
-
y =
-
quadl(@myfun1,1,2)
-
end
x的函式，要能接受向量形式的x輸入，所以構造這個函式的時候考慮到x是向量的情況。
利用arrayfun函式（7.X後的版本都有這個函式，不瞭解這個函式的朋友可以檢視幫助文件，或者百度搜索arrayfun）可以將IntDemo函式精簡成一句程式碼：
-
quadl(@(x)
-
arrayfun(@(xx) quadl(@(y)
- xx*y,sin(xx),cos(xx)),x),1,2)
首先
-
@(x)
-
arrayfun(@(xx) quadl(@(y)
- xx*y,sin(xx),cos(xx)),x)
-
arrayfun(@(xx)
- quadl(@(y) xx*y,sin(xx),cos(xx)),x)
而
-
@(xx) quadl(@(y)
- xx*y,sin(xx),cos(xx))
而
-
arrayfun(@(xx)
- quadl(@(y) xx*y,sin(xx),cos(xx)),x)
-
@(xx) quadl(@(y)
- xx*y,sin(xx),cos(xx))
有了這個模板，我們可以方便求其他一般積分割槽域（上下限是函式）形式的二重積分，例如被積函式
f
= @(x,y)
exp(sin(x))*ln(y)，y從5*x積分到x^2，x從10積分到20。
用quad2d,上面介紹的方法，還有dblquad幫助文件裡給的延拓函式的方法
-
tic,y1
-
= quad2d(@(x,y) exp(sin(x)).*log(y),10,20,@(x)5*x,@(x)x.^2),toc
-
tic,y2 =
-
quadl(@(x) arrayfun(@(x) quadl(@(y)
-
exp(sin(x)).*log(y),5*x,x.^2),x),10,20),toc
-
tic,y3 = dblquad(@(x,y)
-
exp(sin(x)).*log(y).*(y>=5*x & y<=x.^2),10,20,50,400),toc
-
y1
-
=
-
9.368671342614414e+003
-
Elapsed time is 0.021152 seconds.
-
y2
-
=
-
9.368671342161189e+003
-
Elapsed time is 0.276614 seconds.
-
y3
-
=
-
9.368671498376889e+003
-
Elapsed time is 1.674544
-
seconds.