矢量化 - 总和和贝塞尔函数

Question

任何人都可以帮助矢量化此MATLAB代码吗？具体问题是带有向量输入的总和和贝塞尔函数。谢谢！

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

Answer 1

您可以尝试矢量化此代码，这可能在某些bsxfun左右的情况下进行，但是很难理解代码，这是一个问题，这是一个问题，它是否会更快地运行任何速度，由于您的代码已经在内部循环中使用了向量数学（即使您的向量只有长度为3）。由此产生的代码很难阅读，因此您或您的同事在两年内看到它时会做什么。

在浪费矢量化时间之前，了解循环不变的代码运动要易于应用于您的代码，这要重要得多。一些观察：

• 您不使用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的计算需要更多的时间，则可以通过ii和jj交换for-loops，但这似乎并非如此。

结果代码看起来像这样（未经测试）：

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

Answer 2

初始化 -

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。>。这些分别用于矢量化分别计算为术语1,2和3。
3. 这三个术语似乎为1倍大小。因此，我们的目标成为计算这三个术语而无需循环。现在，两个循环每次运行n次，因此为我们提供了NXN的总循环计数。因此，与这些术语在嵌套循环内相比，我们必须在每个术语中具有NXN倍。
。
4. 这基本上是在此处完成的矢量化的本质，因为这三个术语由" term1"，" enter2"one_answers" frc"本身表示。

Answer 3

为了给出一个独立的答案，我将复制原始初始化

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)