0

我正在尝试使用 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, ...
4

2 回答 2

2

你需要一个元素操作符

z=(x0-z0*tan(theta).*cos(phi)).*(y0-z0*tan(theta).*sin(phi))*(z0/2)^4;
                              ^
                              |
于 2013-05-19T14:46:55.307 回答
1

As the help for dblquad says, the input x is a vector and y is a scalar and the output z is also vector. Thus anything in your integrand function that is a function of x will be a vector (e.g., rl) and you'll need to be careful to use element-wise operators where appropriate. You did not do this for the very last line.

Also, consider passing your initial value parameters via function handle rather than duplicating them inside the integrand function:

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;
于 2013-05-19T19:20:05.373 回答