Vectorization of Integral (or quad) to avoid employing a double loop

I wish to vectorize a numeric integration to avoid using a double loop. Here is an example of the problem.

   %%%%%%%%%%%%%%%%%%%%%%%%
% The integral output is actually a function of (x,y) and I am integrating over 
% the variable v, with (x,y) defined using meshgrid.

    [x,y]=meshgrid(-2:0.1:2,-2:0.1,2);
    func=quad( @(v)besselj( 0,v*sqrt(x.^2+y.^2) ),0,1 );
% I would then like to plot the output func as a surface

    mesh(x,y,func);  
%When run, it gives the following error:
??? Error using ==> mtimes
Inner matrix dimensions must agree.

%%%%%%
% I would like to avoid doing this:

    x=linspace(-2,2,21);
    y=linspace(-2,2,21);
    for a=1:1:length(x)
        for b=1:1:length(y)
            func(a,b)=quad(@(v)besselj(0,sqrt(a^2+b^2)*v),0,1);
        end
    end
    mesh(x,y,func);
%Which runs extremely slow.

%%%%%%%%%%%%%%%%%%%%%%%%

When I run the top bit of code I get the following error

Any suggestions?

Answer 1

[Andrew Newell]

You could eliminate one loop by using quadv:

x = -2:0.1:2; y = x;
[X,Y]=meshgrid(-2:0.1:2,-2:0.1,2);

func = zeros(size(X));
for i=1:length(y)
  f = @(v) besselj(0,v*sqrt(x.^2+y(i)^2));
  func(i,:) = quadv(f,0,1);
end
mesh(x,y,func)

This took 0.3 seconds on a Mac Pro.

EDIT: Note that your loops involving a and b give integrals for x from 1 to 21, not x between -2 and 2. I am assuming it is the latter you are really after.

EDIT: Building on Matt's post, which is a superior method but for the wrong problem, here is a very compact solution for your problem:

x = -2:0.1:2; y = x;
[x,y]=meshgrid(x,y);

func = quadv(@(v) besselj(0,v*sqrt(x.^2+y.^2)),0,1);
mesh(x,y,func)

It is now down to 0.08 seconds!

[Matt Fig]

I was wondering about that. Notice the wording of my answer...

Answer 2

[Matt Fig]

This produces the same surface as your double FOR loop:

N = 21;  % The grid

tic  % Your original
x=linspace(-2,2,N);
y=linspace(-2,2,N);
for a=1:1:length(x)
    for b=1:1:length(y)
        func(a,b)=quad(@(v)besselj(0,sqrt(a^2+b^2)*v),0,1);
    end
end
toc

subplot(1,2,1)
mesh(x,y,func)  % Plot your func

tic  % much faster
[x2,y2] = meshgrid(1:N,1:N);
func2 = quadv( @(v)besselj( 0,v.*sqrt(x2.^2+y2.^2) ),0,1 );
toc

subplot(1,2,2)
mesh(linspace(-2,2,N),linspace(-2,2,N),func2);  % Plot new func.

Answer 3

[Walter Roberson]

See the documentation for quad:

"The function y = fun(x) should accept a vector argument x and return a vector result y, the integrand evaluated at each element of x."

Your quad function does not do that: it is written to expect v to be a scalar, and it tries to return an entire array of results.

You will not be able to use quad() as the integrator for the approach you were hoping to use.

You might wish to consider an alternative approach:

int(BesselJ(0,r*v),v=0..1)

is

BesselJ(0, r) - (1/2)*Pi*BesselJ(0, r)*StruveH(1, r) + (1/2)*Pi*BesselJ(1, r)*StruveH(0, r)

I see from the documentation that besselj does accept a matrix as its second argument. Unfortunately, I do not seem to locate a Matlab StruveH . If you have the symbolic toolbox, you could try running an integral similar to what I show above and see what pops up.

The only conversion I have found so far that eliminates the StuveH is via hypergeom:

hypergeom([], [1], -(1/4)*r^2) -
(1/3)*hypergeom([], [1], -(1/4)*r^2) *
r^2*hypergeom([1], [3/2, 5/2], -(1/4)*r^2) +
(1/2)*r^2*hypergeom([], [2], -(1/4)*r^2) *
hypergeom([1], [3/2, 3/2], -(1/4)*r^2)

Answer 4

[Walter Roberson]

Stealing from those who have posted before me ;-)

x = (-2:0.1:2).^2;
r = sqrt(bsxfun(@plus, x.', x));
func = quadv(@(v) besselj(0, v.*r), 0, 1);
mesh(x,x,func);

This is slightly faster than Andrew's code.