0

这是使用 biot-savart-law 计算磁场的代码。我希望得到一些优化此代码的提示。遗憾的是我使用德语:(我再也不会这样做了。:)

tic
clear all; clc; clf
skalierungsfaktor = 10^-6; % vom m-Bereich zum mm-Bereich wg. T = Vs / m^2
I = 1; % in A
konstante = 10^-10; % mu0 / (4 pi) in Vs / (A mm) % Faktor 10^3 kleiner wg. mm statt m
matrix = 30;
B = zeros(200,200,200,5);
zaehler = 0; % für Zeitmessung
xmax = matrix;
ymax = matrix;
zmax = 1;
radius = 9; % Spulenradius in mm
genauigkeit = 5; % 1 = 6.28 Elemente pro Kreis; 2 = 12.56 Elemente pro Kreis; 4 bis 5 scheint gut zu sein
windungen = 10;
leiterelemente = ceil(2 * 3.14152 * genauigkeit * windungen) % 2 * Pi * genauigkeit für eine Umrundung
leiter = repmat(leiterelemente+1,3);
windungen = leiterelemente / genauigkeit / 2 / 3.1415927;
spulenlaenge = 20; % Spulenlaenge in mm
steigung = spulenlaenge / windungen
for i = 1:leiterelemente+1;
    leiter(i,1) = i * steigung / (genauigkeit * 2 * 3.1415927) + matrix/2 - spulenlaenge/2;  % x-Ausrichtung
    leiter(i,2) = radius * cos(i/genauigkeit) + matrix/2;  % y-Ausrichtung
    leiter(i,3) = radius * sin(i/genauigkeit);   % z-Ausrichtung
end
for x = 1:xmax
zaehler = zaehler + 1; % für Zeitmessung
hhh = waitbar(0,num2str(zaehler*100/matrix)); % Wartebalken
waitbar(zaehler/matrix) % Wartebalken
for y = 1:ymax % wenn streamslice nicht genutzt wird, nur einen y-Wert berechnen
    for z = 1:zmax
        for i = 1:leiterelemente
            dl(1) = leiter(i+1,1)-leiter(i,1);
            dl(2) = leiter(i+1,2)-leiter(i,2);
            dl(3) = leiter(i+1,3)-leiter(i,3);
            vecs = [(leiter(i,1)+leiter(i+1,1))/2, ...
                (leiter(i,2)+leiter(i+1,2))/2, ...
                (leiter(i,3)+leiter(i+1,3))/2];
            vecr = [x y z];
            vecrminusvecs = vecr - vecs;
            einheitsvecr = vecrminusvecs./norm(vecrminusvecs); % ok
            r = sqrt(vecrminusvecs(1).^2 + vecrminusvecs(2).^2 + vecrminusvecs(3).^2); % ok
            vektorprodukt = [dl(2).*einheitsvecr(3) - dl(3).*einheitsvecr(2), ...
                dl(3).*einheitsvecr(1) - dl(1).*einheitsvecr(3), ...
                dl(1).*einheitsvecr(2) - dl(2).*einheitsvecr(1)];
            dB = konstante * I * vektorprodukt / (r.^2);
            dB = dB / skalierungsfaktor; % nur hier wird der Wert verändert bzw. skaliert
            B(x,y,z,1) = B(x,y,z,1) + dB(1);
            B(x,y,z,2) = B(x,y,z,2) + dB(2);
            B(x,y,z,3) = B(x,y,z,3) + dB(3);
            B(x,y,z,4) = B(x,y,z,4) + sqrt(dB(1).^2 + dB(2).^2 + dB(3).^2);
        end;
    end;
end;
close(hhh) 
end;
toc
n = 1:leiterelemente;
Lx = leiter(n,1);
Ly = leiter(n,2);
Lz = leiter(n,3);
%subplot(2,1,2), 
line(Lx,Ly,Lz,'Color','k','LineWidth',2); 
hold on
view(15,30);            % view(0,0) = Blickwinkel, 2D-Perspektive
grid on                 % Gitter anzeigen
xlim([0 matrix])
ylim([0 matrix])
zlim([0 5])
xlabel('x-Achse');
ylabel('y-Achse');
zlabel('z-Achse');
daspect([1 1 1])
[X,Y]=meshgrid(1:matrix);
U=(B(1:matrix,1:matrix,z,1))'; 
V=(B(1:matrix,1:matrix,z,2))';
streamslice(X,Y,U,V) % quiver, streamslice
4

4 回答 4

5

不一定是速度优化，但有一些注意事项：

  • 使用pi，而不是 3.14159 或 3.1415927
  • 您通常可以对循环进行矢量化：

非向量化：

for i = 1:leiterelemente+1;
    leiter(i,1) = i * steigung / (genauigkeit * 2 * 3.1415297) + matrix/2 - spulenlaenge/2;  % x-Ausrichtung
    leiter(i,2) = radius * cos(i/genauigkeit) + matrix/2;  % y-Ausrichtung
    leiter(i,3) = radius * sin(i/genauigkeit);   % z-Ausrichtung
end

矢量化：

ii = 1:leiterelemente+1;
leiter(ii,1) = ii * steigung / (genauigkeit * 2 * 3.1415297) + matrix/2 - spulenlaenge/2;  % x-Ausrichtung
leiter(ii,2) = radius * cos(ii/genauigkeit) + matrix/2;  % y-Ausrichtung
leiter(ii,3) = radius * sin(ii/genauigkeit);   % z-Ausrichtung

大多数 matlab 函数将向量/矩阵作为参数，包括cos()、sin()、exp()、log()等。

对于少量元素（比如<几百个），矢量化可能不值得。

向量幅度：而不是sqrt(dB(1).^2 + dB(2).^2 + dB(3).^2)使用norm(dB)（请注意，范数不会以逐行方式对矩阵进行操作，而是在整体上）尽管这不会节省太多

        B(x,y,z,1) = B(x,y,z,1) + dB(1);
        B(x,y,z,2) = B(x,y,z,2) + dB(2);
        B(x,y,z,3) = B(x,y,z,3) + dB(3);

考虑改为

B(x,y,z,1:3) = B(x,y,z,1:3) + dB(1:3)；

为什么在稍后对其进行平方时使用平方根来计算 r？

r = sqrt(vecrminusvecs(1).^2 + vecrminusvecs(2).^2 + vecrminusvecs(3).^2);

改成

r2 = sum(vecrminusvecs.^2);

并使用r2代替r.^2

我的猜测是，您可能可以通过使用一些向量代数来简化从“vecrminusvecs = ...”到“db = konstante ...”的计算；你正在做一些看起来不太必要的重新缩放，或者至少可以针对速度进行优化。

编辑：我现在怀疑“规范”；sqrt(sum(x.^2,2))在每一行上运行，可能比它更快，norm()但如果你想使用最快的方法，你应该测量它。

于 2010-01-25T19:24:00.617 回答
2

你想优化什么？速度？优化的第一个技巧：

测量它。

所以，首先测量你大部分时间花在哪里（可能在某个循环中）。放置计时器，以便您有良好的感觉。然后看看你能不能做点什么。只要你不知道你需要优化什么，你就无法优化。

也就是说，Matlab 的一般规则是不惜一切代价尝试避免循环（尤其是嵌套循环），而是弄清楚如何将计算呈现为矩阵运算。它们快速且经过优化。

于 2010-01-25T17:24:32.360 回答
1

我在 MATLAB 分析器中运行了这段代码，有几件事很突出：

  1. 您每次都在 x = 1:xmax 循环周围创建和销毁等待栏 - 您应该始终保持等待栏打开。这在我的机器上花费了相当多的时间。

  2. 您确实需要至少对三个内部循环进行矢量化。例如，您对 "dl" 的计算可以由对 (something like - untested) 的单个调用替换dl = diff( leiter, 1, 1 )。同样，vecs = (leiter(1:N-1,:) + leiter(2:N,:))/2. 该循环中的其余表达式需要更多的工作来梳理。

于 2010-01-26T08:21:28.410 回答
0

这只是一个 20ppm 的错误，但 pi ~= 3.1415927 != 3.1415297。（你有一个错字，切换了 2 和 9）。除了方便使用内置常量之外的另一个原因。

于 2010-03-11T04:08:02.933 回答