Numerical integration with array limits

I'm trying to work with integrals that are functions of one of their limits:

For example,

phi = @(x) quad(@(L) besseli(1, (1+L)/(1-L)), 0, x);

What I'm trying to do is evaluate phi over an array of values, like:

phi([1,2,3,4]); %ERROR
quad(@(L) besseli(1, (1+L)/(1-L)), 0, [1,2,3,4]); %ERROR

but these return errors. I could do this in a for loop, like:

nums=[1,2,3,4];
for(k=1:4)
  phi_eval = phi(nums(k));
end

but I was wondering if there was a better way to do things. Is there a no-for-loops way of doing this?

Answer 1

[Shashank Prasanna]

Christopher, I can't run the loop as well. But phi([1,2,3,4]) will certainly not work because the vector is being passed to quad directly as limits which is wrong syntax for quad.

You can try the following:

arrayfun(phi,[1,2,3,4])

http://www.mathworks.com/help/matlab/ref/arrayfun.html

[Christopher Kodama]

Great, thanks! If I had known about arrayfun, I probably would've done that from the beginning.

I also found another solution that involved u-substitution in the integral (which allowed the limits to be from 0 to 1 instead of 0 to L) and the quadv function:

nums = [1,2,3,4];
phi = quadv(@(L) besseli(1, (1+L.*nums)./(1-L.*nums)).*nums, 0, 1);

[Mike Hosea]

Note that both quad and quadv are deprecated. Use integral instead. If you want to integrate a vector-valued problem, use

integral(...,'ArrayValued',true)

[Wenjuan]

What if both limits are vectors as well? I don't think integral, quad or quadv can deal with this, but how to use arrayfun in this? Thanks.

[Leo Simon]

I too would like to specify vector valued integration limits. If anybody from mathworks is listening could you please respond? This should be really easy to implement I imagine, and hopefully will be in the next release?

[Mike Hosea]

Please go to my profile. Where it says "email", click on "contact Mike Hosea" and tell me about your use cases for array-valued limits. If we're talking generic array limits, where there is no a priori relationship between the different elements of the limits, then no gain in efficiency can be had over writing a loop. E.g.

Q = zeros(size(a));
for k = 1:numel(Q)
    Q(k) = integral(f,a(k),b(k));
end

For scalar-valued integrations, that can also be accomplished efficiently with an application of arrayfun, e.g.

Qarray = @(a,b)arrayfun(@(ak,bk)integral(f,ak,bk),a,b);
Q = Qarray(a,b);

However, if we are talking about table-building, where the limits represent a grid, then efficiency improvements are possible. The latter use case might be accomplished by some other means, however, such as providing scalar limits and a list of output points for the integrals over partial regions.

[Ashish Bhatt]

it works with both array limits as well:

phi = matlabFunction(int(x^2,x,a,b));
arrayfun(phi,vec1,vec2)