Numerical Integration in MatLab

Question

I am trying to compute an integral using Matlab's dblquad. To do it I first wrote a script function. My code is

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;

To compute the numerical integral I type in command window

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

I get an error saying that

??? 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, ...

What am I doing wrong with that?

EDIT

Replacing

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

with

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

I get an new error

??? 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

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;

Answer 2

You need an element-wise operator

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