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

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

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)

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);

y = f(x, varargin{:});

Q(i) = quadf(intfcn, xmin, xmax, tol, trace, y(i), varargin{:});

y = f(x, varargin{:});

Q = quadf(@innerintegral, ymin, ymax, tol, trace, intfcn, ...

Error in I2 (line 8)


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)

-
If I understand correctly, looks like you didn't set a value for sigma when calling Hsigma. You need to change it to something like out1=H(r0,r1,1) (for sigma=1). When you think about, when you integrate a function of three variables over two, you need to set a constant value for the third one. – Gil Feb 17 '13 at 19:23

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.

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

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);

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.

-
So if you have a nested function g inside a function f and a variable is defined in f it is automatically accesable in g ? – k1next Feb 17 '13 at 20:29
@macydanim Yes, unless you specify it as a parameter to the inner function. See this section from the Mathworks MATLAB page on nested functions. What you want to do is similar to their parabola example. – Daryl Feb 17 '13 at 20:38
@StefanSmith dblquad expects that the function can evaluate in a specific format. In the documentation, it specifies that the function should be able to evaluate with the first input (I.e. r0) as a vector and the second input as a scalar. I suspect this is now your problem. To overcome this, you need to rewrite H.m to accept a vector input for r0. – Daryl Feb 18 '13 at 2:17
@StefanSmith Please see the edit(s) for additional (hopefully helpful) information which I believe will help solve your problem. – Daryl Feb 18 '13 at 3:48
@Daryl : I am still having some trouble. I looked at the documentation link that you give, and it still wasn't clear to me. I found a page (people.whitman.edu/~hundledr/courses/M467/MatlabQuad.pdf) that seems to explain it, but I would like confirmation. Please tell me if this is correct: in the integrand function, Matlab expects the first argument to be a vector and the rest to be scalars. The integrand function must return a vector of the same length as the first input. You can loop through the components of the first input, or use Matlab's vectorising abilities. – Stefan Smith Feb 22 '13 at 0:54

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;



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.

-