Vectorization of double for loop including sine of two variables

Question

I need to numerically evaluate some integrals which are all of the form shown in this image:

These integrals are the matrix elements of a N x N matrix, so I need to evaluate them for all possible combinations of n and m in the range of 1 to N. The integrals are symmetric in n and m which I have implemented in my current nested for loop approach:

function [V] = coulomb3(N, l, R, R0, c, x)
r1 = 0.01:x:R;
r2 = R:x:R0;
r = [r1 r2];
rl1 = r1.^(2*l);
rl2 = r2.^(2*l);
sines = zeros(N, length(r));
V = zeros(N, N);
for i = 1:N;
    sines(i, :) = sin(i*pi*r/R0);
end

x1 = length(r1);
x2 = length(r);
for nn = 1:N
    for mm = 1:nn
        f1 = (1/6)*rl1.*r1.^2.*sines(nn, 1:x1).*sines(mm, 1:x1);
        f2 = ((R^2/2)*rl2 - (R^3/3)*rl2.*r2.^(-1)).*sines(nn, x1+1:x2).*sines(mm, x1+1:x2);
        value = 4*pi*c*x*trapz([f1 f2]);
        V(nn, mm) = value;
        V(mm, nn) = value;
    end
end

I figured that calling sin(x) in the loop was a bad idea, so I calculate all the needed values and store them. To evaluate the integrals I used trapz, but as the first and the second/third integrals have different ranges the function values need to be calculated separately and then combined.

I've tried a couple different ways of vectorization but the only one that gives the correct results takes much longer than the above loop (used gmultiply but the arrays created are enourmous). I've also made an analytical solution (which is possible assuming m and n are integers and R0 > R > 0) but these solutions involve a cosine integral (cosint in MATLAB) function which is extremely slow for large N.

I'm not sure the entire thing can be vectorized without creating very large arrays, but the inner loop at least should be possible. Any ideas would be be greatly appreciated!

The inputs I use currently are:

R0 = 1000;
R = 8.4691;
c = 0.393*10^(-2);
x = 0.01;
l = 0 # Can reasonably be 0-6;
N = 20; # Increasing the value will give the same results, 
# but I would like to be able to do at least N = 600;

Using these values

V(1, 1:3) = 873,379900963549    -5,80688363271849   -3,38139152472590

Although the diagonal values never converge with increasing R0 so they are less interesting.

Answer 1

You will lose the gain from the symmetricity of the problem with my approach, but this means a factor of 2 loss. Odds are that you'll still benefit in the end.

The idea is to use multidimensional arrays, making use of trapz supporting these inputs. I'll demonstrate the first term in your figure, as the two others should be done similarly, and the point is the technique:

r1 = 0.01:x:R;
r2 = R:x:R0;
r = [r1 r2].';
rl1 = r1.'.^(2*l);
rl2 = r2.'.^(2*l);
sines = zeros(length(r),N);       %// CHANGED!!
%// V = zeros(N, N);  not needed now, see later
%// you can define sines in a vectorized way as well:
sines = sin(r*(1:N)*pi/R0);     %//' now size [Nr, N] !

%// note that implicitly r is of size [Nr, 1, 1]
%// and sines is of size [Nr, N, 1]
sines2mat = permute(sines,[1, 3, 2]);  %// size [Nr, 1, N]

%// the first term in V: perform integral along first dimension
%//V1 = 1/6*squeeze(trapz(bsxfun(@times,bsxfun(@times,r.^(2*l+2),sines),sines2mat),1))*x;  %// 4*pi*c prefactor might be physics, not math
V1 = 1/6*permute(trapz(bsxfun(@times,bsxfun(@times,r.^(2*l+2),sines),sines2mat),1),[2,3,1])*x;  %// 4*pi*c prefactor might be physics, not math

The key point is that bsxfun(@times,r.^(2*l+2),sines) is a matrix of size [Nr,N,1], which is again multiplied by sines2mat using bsxfun, the result is of size [Nr,N,N] and an element (k1,k2,k3) corresponds to an integrand at radial point k1, n=k2 and m=k3. Using trapz() with explicitly the first dimension (which would be default) reduces this to an array of size [1,N,N], which is just what you need after a good squeeze(). Update: as per @Dev-iL's comment you should use permute instead of squeeze to get rid of the leading singleton dimension, as that might be more efficent.

The two other terms can be handled the same way, and of course it might still help if you restructure the integrals based on overlapping and non-overlapping parts.