1

任何人都可以帮助矢量化这个 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
4

3 回答 3

1

您可以尝试对这段代码进行矢量化，这可能有一些可能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
于 2014-03-07T20:10:12.110 回答
0

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

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)];

并生成一些缺失数据（n维的k(3)和rho_s和phi_s）

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)
于 2014-03-07T23:00:15.680 回答
0

初始化 -

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));

关于这种矢量化或代码简化的注意事项 –</p>

  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”本身表示。
于 2014-03-07T23:07:58.883 回答