Векторизация двойного цикла, включая синус два переменных

Question

мне нужно численно оценить некоторые интегралы, которые все форм, показанных на этом изображении:

Этих интегралов являются матричными элементами N x N матрицы, так что мне нужно для их оценки для всех возможных комбинаций n и m в диапазоне от 1 до N. Интегралы симметричны n и m, которые я реализовал в моем текущем вложенной подходе for петли:

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 

Я понял, что вызов sin(x) в петле была плохая идея, так что я вычислить все необходимые значения и хранить их , Чтобы оценить интегралы, я использовал trapz, но поскольку первый и второй/третий интегралы имеют разные диапазоны, значения функций нужно вычислять отдельно и затем объединяться.

Я пробовал пару различных способов векторизации, но единственный, который дает правильные результаты, занимает намного больше времени, чем приведенный выше цикл (используется gmultiply, но созданные массивы являются излишними). Я также сделал аналитическое решение (что можно предположить, что m и n являются целыми числами и R0 > R > 0), но эти решения включают интеграл косинуса (cosint в MATLAB), который чрезвычайно мал для больших N.

Я не уверен, что вся вещь может быть векторизована без создания очень больших массивов, но по крайней мере должна быть внутренняя петля. Любые идеи были бы очень признательны!

входы, которые я использую в настоящее время являются:

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; 

Используя эти значения

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

Хотя диагональные значения никогда не сходятся с увеличением R0 поэтому они менее интересны.

Answer 1

Вы потеряете выигрыш от симметричности проблемы с моим подходом, но это означает коэффициент потери 2. Скорее всего, вы все равно выиграете в конце.

Идея заключается в использовании многомерных массивов с использованием trapz, поддерживающих эти входы. Я продемонстрирую первый член в вашей фигуре, а два других должно быть сделано так же, и точка является метод:

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 

Ключевым моментом является то, что bsxfun(@times,r.^(2*l+2),sines) представляет собой матрицу размером [Nr,N,1], что опять-таки умноженное на sines2mat с использованием bsxfun, результат будет иметь размер [Nr,N,N], а элемент (k1,k2,k3) соответствует подынтегральному выражению в радиальной точке k1, n=k2 и m=k3. Используя trapz() с явно первым измерением (которое было бы по умолчанию), это уменьшит это до массива размера [1,N,N], который вам нужен только после хорошего squeeze(). Обновление: согласно @Dev-iL's comment, вы должны использовать permute вместо squeeze, чтобы избавиться от основного размера синглтона, поскольку это может быть более эффективным.

Два других терминала могут обрабатываться одинаково, и, конечно, это может помочь, если вы реструктурируете интегралы на основе перекрывающихся и неперекрывающихся частей.