Inconsistency using integral versus quad?

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?

Answer 1

[Mike Hosea]

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.

[Mike Hosea]

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.

Answer 2

[Sam Butler]

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.