Векторизация — сумма и функция Бесселя

Инициализация -

N = 3;
rho_g = linspace(1e-3,1,N);
phi_g = linspace(0,2*pi,N);

n = 1:3;
tau = [1 2.*ones(1,length(n)-1)];

Форма вложенных циклов (скопировано из вашего кода и показано здесь только для сравнения) -

for ii = 1:length(rho_g)
    for jj = 1:length(phi_g)
        % Coordinates
        rho_o = rho_g(ii);
        phi_o = phi_g(jj);
        % factors
        fc = cos(n.*(phi_o-phi_s));
        fs = sin(n.*(phi_o-phi_s));

        Ez_t(ii,jj) = sum(tau.*besselj(n,k(3)*rho_s).*besselh(n,2,k(3)*rho_o).*fc);
    end
end

Векторное решение -

%%// Term - 1
term1 = repmat(tau.*besselj(n,k(3)*rho_s),[N*N 1]);

%%// Term - 2
[n1,rho_g1] = meshgrid(n,rho_g);
term2_intm = besselh(n1,2,k(3)*rho_g1);
term2 = transpose(reshape(repmat(transpose(term2_intm),[N 1]),N,N*N));

%%// Term -3
angle1 = repmat(bsxfun(@times,bsxfun(@minus,phi_g,phi_s')',n),[N 1]);
fc = cos(angle1);

%%// Output
Ez_t = sum(term1.*term2.*fc,2);
Ez_t = transpose(reshape(Ez_t,N,N));

Замечания об этой векторизации или упрощении кода:

1. «fs» не меняет вывод скрипта, Ez_t, поэтому его можно пока удалить.
2. Вывод выглядит как «Ez_t», что требует трех основных терминов в коде: tau.*besselj(n,k(3)*rho_s), besselh(n,2, k(3)*rho_o) и fc. Они рассчитываются отдельно для векторизации как условия 1, 2 и 3 соответственно.
3. Все эти три термина имеют размеры 1xN. Таким образом, наша цель состоит в том, чтобы вычислить эти три члена без циклов. Теперь два цикла выполняются по N раз каждый, что дает нам общее количество циклов NxN. Таким образом, мы должны иметь в NxN раз больше данных в каждом таком терме по сравнению с тем, когда эти термы находились внутри вложенных циклов.
4. Это в основном суть векторизации, выполненной здесь, поскольку три термина представлены как «term1», «term2» и «fc».

Вы можете попытаться векторизовать этот код, что может быть возможно с некоторыми bsxfun или около того, но будет трудно понять код, и вопрос в том, будет ли он работать быстрее, поскольку ваш код уже использует векторную математику во внутреннем цикле. (хотя ваши векторы имеют только длину 3). Полученный код станет очень трудным для чтения, поэтому вы или ваш коллега не будете знать, что он делает, когда взглянете на него через 2 года.

Прежде чем тратить время на векторизацию, гораздо важнее узнать о движении кода, инвариантном к циклу. , который легко применить к вашему коду. Некоторые наблюдения:

• вы не используете fs, поэтому удалите это.
• термин tau.*besselj(n,k(3)*rho_s) не зависит ни от одной из ваших переменных цикла ii и jj, поэтому он постоянен. Вычислите его один раз перед циклом.
• вам, вероятно, следует предварительно выделить матрицу Ez_t.
• единственные термины, которые изменяются во время цикла, это fc, который зависит от jj, и besselh(n,2,k(3)*rho_o), который зависит от ii. Я предполагаю, что последнее требует гораздо больше времени для вычисления, поэтому лучше не вычислять это N*N раз во внутреннем цикле, а только N раз во внешнем цикле. Если бы вычисление на основе jj заняло больше времени, вы могли бы поменять местами циклы for для ii и jj, но здесь это не так.

Код результата будет выглядеть примерно так (не проверено):

N = 3;
rho_g = linspace(1e-3,1,N);
phi_g = linspace(0,2*pi,N);

n = 1:3;
tau = [1 2.*ones(1,length(n)-1)];

% constant part, does not depend on ii and jj, so calculate only once!
temp1 = tau.*besselj(n,k(3)*rho_s); 

Ez_t = nan(length(rho_g), length(phi_g)); % preallocate space
for ii = 1:length(rho_g)
    % calculate stuff that depends on ii only
    rho_o = rho_g(ii);
    temp2 = besselh(n,2,k(3)*rho_o);

    for jj = 1:length(phi_g)
        phi_o = phi_g(jj);
        fc = cos(n.*(phi_o-phi_s));
        Ez_t(ii,jj) = sum(temp1.*temp2.*fc);
    end
end

Чтобы дать самостоятельный ответ, я скопирую исходную инициализацию

N = 3;
rho_g = linspace(1e-3,1,N);
phi_g = linspace(0,2*pi,N);

n = 1:3;
tau = [1 2.*ones(1,length(n)-1)];

и сгенерировать некоторые недостающие данные (k(3) и rho_s и phi_s в измерении n)

rho_s = rand(size(n));
phi_s = rand(size(n));
k(3) = rand(1);

то вы можете вычислить тот же Ez_t с многомерными массивами:

[RHO_G, PHI_G, N] = meshgrid(rho_g, phi_g, n);
[~, ~, TAU] = meshgrid(rho_g, phi_g, tau);
[~, ~, RHO_S] = meshgrid(rho_g, phi_g, rho_s);
[~, ~, PHI_S] = meshgrid(rho_g, phi_g, phi_s);

FC = cos(N.*(PHI_G - PHI_S));
FS = sin(N.*(PHI_G - PHI_S)); % not used

EZ_T = sum(TAU.*besselj(N, k(3)*RHO_S).*besselh(N, 2, k(3)*RHO_G).*FC, 3).';

После этого вы можете проверить, что обе матрицы одинаковы.

norm(Ez_t - EZ_T)