Inconsistency using integral versus quad?

5 views (last 30 days)
Sam Butler
Sam Butler on 25 Aug 2014
Edited: Mike Hosea on 5 Sep 2014
I am encountering a strange problem using the built-in "integral" function in MatLab R2013b, where it is not giving me the proper answer. I am wondering if it is a usage difference between integral and quad, which seem at a glance to be very similar (other than the suggestion to use "integral" instead, which is where the bug resides).
I have a function which is a normalized Gaussian in theta-space, so that if I integrate it over all angular coordinates in a hemisphere, it should integrate to 1...that is, 2 * pi * integral(f(t)*cos(t)*sin(t), t=0 to pi/2) = 1 since it is azimuthally symmetric. (The cos(t) is due to power projection, sin(t) comes from spherical coordinate integration. The 2*pi is the integration over phi due to symmetry.)
I declare an inline function: intfunc = @(t) f(t) .* cos(t) .* sin(t)
Then I tried the following two snippets of code: intval = 2*pi*integral(intfunc, 0, pi/2)
The result that I get varies depending upon the choice of Gaussian width, but it does not integrate to 1.
If instead I do the following: intval = 2*pi*quad(intfunc, 0, pi/2)
Then the result integrates to 1. So, it appears one of the functions isn't working the way I expect.
To resolve this dilemma, I went to Mathematica, integrated the function symbolically, and obtained a value of 1 for any value of the Gaussian width (as long as it is positive, of course). This suggests to me that quad works properly and integral does not.
Am I somehow misusing integral(), or is it known to be buggy?

Accepted Answer

Mike Hosea
Mike Hosea on 2 Sep 2014
Edited: Mike Hosea on 2 Sep 2014
No, INTEGRAL is extremely reliable, much better than QUAD. Please provide me with a complete, running example that demonstrates the problem you are having (the smaller and more meaningless the better :) ), and I'll be happy to sort out what's going on.
The most common problems occur when the integrand has not been defined to operate on vector inputs. Per the documentation, your integrand function must be able to accept an array input (of arbitrary size, though I think it is always a vector), evaluate the integrand for each element of the input array, and return the corresponding array of outputs. The most common errors are to use a * instead of .*, ^ instead of .^, and / instead of ./. If it is difficult to "vectorize" an integrand in this manner, then you have a couple of options. One is to integrate @(x)arrayfun(f,x) instead of just f. The arrayfun "wrapper" just takes an array and evaluates the function f "element-wise". Another way is to use the 'ArrayValued',true option to INTEGRAL. This is intended for integrating vector-valued integrands, i.e. integrands that return a vector (or some other array shape) when you give them a scalar. However, a scalar return is a perfectly valid "array". Consequently, 'ArrayValued',true has the effect of turning off the vectorized-integrand requirement.
Of course I see that your intfun definition above is vectorized, but what I can't see is whether f(t) is.
  1 Comment
Mike Hosea
Mike Hosea on 5 Sep 2014
Thanks. Vectorization was, indeed, the issue. I don't know what "vec_half" is all about there, but that's your problem. Regardless of how many x-values INTEGRAL sent in, you were only ever looking at the first one. This should work better
function [val] = Dist_Gaussian(theta, params)
sig = params(1);
assert(sig > 0);
prefix = (2.*pi.*sig.^2.*cos(theta).^4);
val = (1 ./ prefix) .* exp(-0.5 .* (tan(theta)./sig).^2 );
end
As I said before, if you have a problem that is more complicated and more difficult to fix, you can integrate @(x)arrayfun(f,x) instead or tack on the 'ArrayValued',true to the end of the INTEGRAL call. The fastest performance, of course, will be had if f can be vectorized naturally.

Sign in to comment.

More Answers (1)

Sam Butler
Sam Butler on 3 Sep 2014
Here is f(x); it is actually a two-parameter code:
function [val] = Dist_Gaussian(vec_half, params) sig = params(1); assert(sig > 0); theta = vec_half(1); prefix = (2.*pi.*sig.^2.*cos(theta).^4); val = (1 ./ prefix) .* exp(-0.5 .* (tan(theta)./sig).^2 ); end
To integrate it, I use:
p = [0.1]; f = @(x) Dist_Gaussian(x, p) .* cos(x) .* sin(x);
Then I integrate this as per above.
This will give the right answer using quad, but not using integral.

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!