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;