MatLab中的數值積分

Question

我想用Matlab的dblquad來計算積分。要做到這一點，我首先寫了一個腳本函數。我的代碼是

function z = IntRect(x, y) 
%The variables 
z0=20; 
x0=15; 
y0=20; 
rl=sqrt(x.^2+y.^2+(z0/2)^2); 
theta=acos(z0./(2*rl)); 
phi=atan(y./x); 
%The function 
z=(x0-z0*tan(theta).*cos(phi))*(y0-z0*tan(theta).*sin(phi))*(z0/2)^4; 
z=z/rl.^3; 

要計算數值積分我在命令窗口

z0=20;x0=15;y0=20; 
Q = dblquad(@IntRect,0,x0/2,0,y0/2,1.e-6); 

型我得到一個錯誤，說

??? Error using ==> mtimes 
Inner matrix dimensions must agree. 

Error in ==> IntRect at 8 
z=(x0-z0*tan(theta).*cos(phi))*(y0-z0*tan(theta).*sin(phi))*(z0/2)^4; 

Error in ==> quad at 77 
y = f(x, varargin{:}); 

Error in ==> dblquad>innerintegral at 84 
    Q(i) = quadf(intfcn, xmin, xmax, tol, trace, y(i), varargin{:}); 

Error in ==> quad at 77 
y = f(x, varargin{:}); 

Error in ==> dblquad at 60 
Q = quadf(@innerintegral, ymin, ymax, tol, trace, intfcn, ... 

我在做什麼不好呢？

編輯

更換

z=(x0-z0*tan(theta).*cos(phi))*(y0-z0*tan(theta).*sin(phi))*(z0/2)^4; 

與

z=(x0-z0*tan(theta).*cos(phi)).*(y0-z0*tan(theta).*sin(phi))*(z0/2)^4; 

我得到一個新的錯誤

??? Index exceeds matrix dimensions. 

Error in ==> quad at 85 
if ~isfinite(y(7)) 

Error in ==> dblquad>innerintegral at 84 
    Q(i) = quadf(intfcn, xmin, xmax, tol, trace, y(i), varargin{:}); 

Error in ==> quad at 77 
y = f(x, varargin{:}); 

Error in ==> dblquad at 60 
Q = quadf(@innerintegral, ymin, ymax, tol, trace, intfcn, ... 

Answer 1

由於幫助dblquad說，輸入x是一個載體，y是一個標量和輸出z也向量。因此，您的被積函數中的任何函數是x的函數將是一個向量（例如，rl），您需要在適當的地方使用基於元素的運算符。最後一行你沒有這樣做。

另外，還要考慮通過函數句柄傳遞你的初始值的參數，而不是複製它們被積函數內部：

function dblquadtest 
z0 = 20; x0 = 15; y0 = 20; 
f = @(x,y)IntRect(x,y,x0,y0,z0); 
Q = dblquad(f,0,x0/2,0,y0/2); % 1e-6 is the default tolerance 

function z = IntRect(x, y, x0, y0, z0) 
%The variables 
rl=sqrt(x.^2+y.^2+(z0/2)^2); 
theta=acos(z0./(2*rl)); 
phi=atan(y./x); 
%The function 
z=(x0-z0*tan(theta).*cos(phi)).*(y0-z0*tan(theta).*sin(phi))*(z0/2)^4; 
z=z./rl.^3; 

Answer 2

你需要一個eleme NT-明智的運營商

z=(x0-z0*tan(theta).*cos(phi)).*(y0-z0*tan(theta).*sin(phi))*(z0/2)^4; 
          ^
           |