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

你需要一个元素操作符

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

Answer 2

正如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;