Numerically integrating a function of 3 variables with respect to 2 variables in Matlab

Question

Math people:

I have Googled this question, tried some suggestions, and none of them worked. I browsed the Similar Questions here and didn't find what I wanted. I apologize if this is a duplicate. I posted the question on a Matlab forum, but questions on Math Stack Exchange get answered EXREMELY quickly, so I will try here.

I have a function H defined explicitly by algebraic formulas and the heaviside function in an .m file. The function works. I am trying to write a function that integrates H(r0,r1,sigma) numerically as (r0,r1) ranges over a square. Here is my attempt, using a nested function. I'm sorry, but I don't know how to make this look nice in $\LaTeX$. Part of it is in a gray box, and part of it isn't. Some line breaks are missing.

function out = I2(sigma,Rmax) 
% H is a function of 3 variables defined by an m-file, that m-file works.

function out1 = Hsigma(r0,r1)%nested function
    out1 = H(r0,r1,sigma);
end   

out = dblquad(Hsigma,1,Rmax,1,Rmax);
end

I get the following error (again, some of it is in a gray box:

>> I2(2,2)
Error using I2/Hsigma (line 6)
Not enough input arguments.

Error in I2 (line 10)
  out = dblquad(Hsigma,1,Rmax,1,Rmax);

NOTE: line 6 is "out1 = H(r0,r1,sigma); ".

Certainly people have done this before and someone out there knows how to do this. Any ideas?

UPDATE: It looks like someone edited my question to make it look nicer. Thanks. It tried changing the line "out = dblquad(Hsigma,1,Rmax,1,Rmax);" to "out = dblquad(@Hsigma,1,Rmax,1,Rmax);" and I get the following bizarre error messages, starting with a complaint about my function H, which works only with real numbers and not with vectors or matrices:

>> I2(2,2)
Error using  * 
Inner matrix dimensions must agree.

Error in H (line 5)
out = r0*r0*r0/(r1*r1*r1)*heaviside(r0-r1)*2* ...

Error in I2/Hsigma (line 6)
        out1 = H(r0,r1,sigma);

Error in quad (line 72)
y = f(x, varargin{:});

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

Error in quad (line 72)
y = f(x, varargin{:});

Error in dblquad (line 58)
Q = quadf(@innerintegral, ymin, ymax, tol, trace, intfcn, ...

Error in I2 (line 8)
  out = dblquad(@Hsigma,1,Rmax,1,Rmax);

H.m is a pretty innocuous function file. Here it is. I cut-and-pasted it from the Matlab editor, and the * symbol does not appear in the question some places where it should:

function out= H(r0,r1,sigma)
%integrand for I2.  Contains
%integral over an infinite integral, computed by Maple
%   
out = r0*r0*r0/(r1*r1*r1)*heaviside(r0-r1)*2* ...
    ((exp(-(1/4)*(r0^2+2*r0*r1+r1^2)/sigma^2)*r0 ...
        +exp(-(1/4)*(r0^2+2*r0*r1+r1^2)/sigma^2)*r1 ...
        -pi^(1/2)*erf((1/2)*(r0+r1)/sigma)*sigma ...
        +pi^(1/2)*sigma)/sigma);
end

I apologize again for not knowing how to format this correctly. I am at a loss. Can anyone help?

ANOTHER UPDATE: I tried doing the same thing using the "dummy function" Hdummy(r0,r1,sigma)=r0+r1+sigma and it seems to work. I don't know why it doesn't work with the more complicated function H.

Stefan (STack Exchange FAN)

Answer 1

The line out = dblquad(Hsigma,1,Rmax,1,Rmax); is incorrect. MATLAB thinks you are calling Hsigma here, rather than passing a function handle. I.e. it thinks you want to call Hsigma with no parameters. If you change the line to

out = dblquad(@Hsigma,1,Rmax,1,Rmax);

everything should be sorted.

To use a function as an input, you need to specify the @function, which is MATLAB's way of speceifying a function handle.

The problem is not with sigma. It will be defined in the outer function, and shared to the inner function with MATLAB's automatic nested function variable sharing.

Added:

There are two ways which you can achieve the desired functionality. One way is the method which you have used, and I will include below. Note the difference of the @Hsigma on line 7. This tells MATLAB that you are passing the function handle to dblquad and not evaluating the function.

function out = I2(sigma,Rmax)

    function out1 = Hsigma(r0,r1)
        out1 = H(r0,r1,sigma);
    end

    out = dblquad(@Hsigma,1,Rmax,1,Rmax);

end

The second method is to define an anonymous function in the script. It operates in the same fashion as this does, but the syntax is slightly different. Note that Hsigma2 is a function handle to the anonymous function @(r0,r1) H(r0,r1,sigma). Since Hsigma2 is already a function handle, no @ symbol is required in the call to dblquad.

function out = I2(sigma,Rmax)

    Hsigma2 = @(r0,r1) H(r0,r1,sigma);

    out = dblquad(Hsigma2,1,Rmax,1,Rmax);

end

In both cases, the value used for sigma is obtained from the input to the function I2. Both cases are also handled similarly in MATLAB.

Edit

As per my comments below, dblquad expects the function to be able to handle a vector as the first argument. One way to reqrite H.m to achieve this is below. Note that this is probably not the most efficient way to do this, but it is the easiest given your current work.

function out= H(r0vec,r1,sigma)
%integrand for I2.  Contains
%integral over an infinite integral, computed by Maple
%   

n = length(r0vec); % get the length of the input vector
out = zeros(n,1); % allocate memory for output array

for i = 1:n % loop over the vector array
    r0 = r0vec(i); % get the relevant index out of the array
    out(i) = r0*r0*r0/(r1*r1*r1)*heaviside(r0-r1)*2* ...
               ((exp(-(1/4)*(r0^2+2*r0*r1+r1^2)/sigma^2)*r0 ...
               +exp(-(1/4)*(r0^2+2*r0*r1+r1^2)/sigma^2)*r1 ...
               -pi^(1/2)*erf((1/2)*(r0+r1)/sigma)*sigma ...
               +pi^(1/2)*sigma)/sigma);
end % end of loop i

end

All the changes to this function do is loop over each of the components in the first argument to the function. An alternative (and computationally faster) method is to vectorise each operation involving r0.

Edit again

The function is able to be vectorised relatively trivially in MATLAB.

function out= H(r0,r1,sigma)
%integrand for I2.  Contains
%integral over an infinite integral, computed by Maple

    PC = (r0+r1)/(2*sigma); % precalculate this quantity
    out = heaviside(r0-r1).*(r0./r1).^3*2/sigma.*...
            ((r0+r1).*exp(-1*PC.^2)+sqrt(pi)*sigma*(1-erf(PC)));
end

This differs primarily from the previous form in the use of MATLAB's . operators, which indicate a vector operation, rather than a regular operation. e.g. A.*B multiplies corresponding elements of A and B, whereas A*B performs the matrix multiplication of A and B.

Furthermore, this form is vectorised in both r0 and r1, which allows the use of the MATLAB function quad2d. The calling syntax is similar to dblquad. Further information about quad2d can be found here.

Answer 2

As Gil already mentioned in his comment you don't specify a value for sigma anywhere. You could rewrite your code as

function out1 = Hsigma(r0,r1)%nested function
    sigma = 1;
    out1 = H(r0,r1,sigma);
end

And then call

out = dblquad(Hsigma,1,Rmax,1,Rmax);

As and alternative you could try whether the following code also works with dblquad.

function y = myfunc(x,c)
    y = c*x
end

c=2;

quad(@(x) myfunc(x,c), 0, 1.0)

EDIT: From this page:

dblquad(fun,xmin,xmax,ymin,ymax,tol,method,p1,p2,...)

or with default values

dblquad(fun,xmin,xmax,ymin,ymax,[],[],p1,p2,...)

which is equivalent to

dblquad(fun,xmin,xmax,ymin,ymax,1.e-6,@quad,p1,p2,...)

EDIT2: As pointed out by Daryl it has to be @Hsigma instead of Hsigma.