matlab - Matlab：biot-savarts-lawで計算された磁場

Question

これは、ビオ・サバールの法則を使用して磁場を計算するためのコードです。このコードを最適化するためのヒントをいくつか入手したいと思います。残念ながら私はドイツ語を使用しています:(これは二度と行いません。:)

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

score 5 · Accepted Answer

必ずしも速度の最適化ではありませんが、いくつかの注意事項:

pi3.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()が、最速のアプローチを使用する場合は測定する必要があります。

score 2 · Accepted Answer

何を最適化したいですか？スピード？最適化の最初のヒント：

それを測定します。

したがって、最初に、ほとんどの時間を（おそらくいくつかのループで）費やす場所を測定します。あなたがそれをよく感じるようにタイマーを置いてください。次に、それについて何かできるかどうかを確認します。何を最適化する必要があるかわからない限り、最適化することはできません。

とは言うものの、Matlabの一般的なルールは、ループ（特にネストされたループ）を絶対に避け、代わりに行列演算として計算を提示する方法を理解することです。それらは高速で最適化されています。

score 1 · Accepted Answer

このコードを MATLAB プロファイラーで実行したところ、いくつかの点が際立っていました。

x = 1:xmax ループのたびにウェイトバーを作成および破棄しています。ウェイトバーは常に開いたままにしておく必要があります。これは、私のマシンではかなりの時間を要していました。
少なくとも 3 つの内部ループをベクトル化する必要があります。たとえば、「dl」の計算は、 (something - untested) への単一の呼び出しに置き換えることができますdl = diff( leiter, 1, 1 )。同様に、vecs = (leiter(1:N-1,:) + leiter(2:N,:))/2. そのループ内の残りの式を分解するには、もう少し作業が必要です。

score 0 · Accepted Answer

これはわずか 20ppm のエラーですが、pi ~= 3.1415927 != 3.1415297 です。（2と9を入れ替えたタイプミスがあります）。組み込み定数を使用する便利さ以外のもう 1 つの理由。

matlab - Matlab：biot-savarts-lawで計算された磁場

4 に答える 4

Related

Reference