Matlab: магнитные поля, рассчитанные по закону Био-Саварта

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

Answer 1

не обязательно оптимизация скорости, но некоторые примечания:

• используйте pi, а не 3.14159 или 3.1415927
• обычно вы можете векторизовать свои циклы:

nonvectorized:

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() и т. Д.

Для небольшого количества элементов (скажем, <несколько сотен) может быть не стоит усилий по векторизации. </p>

Величина вектора: вместо 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) + дБ (1: 3);

Почему вы вычисляете r с использованием квадратного корня, когда вы просто возводите его в квадрат позже?

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

Изменить на

r2 = sum(vecrminusvecs.^2);

и используйте r2 вместо r.^2

Полагаю, вы, вероятно, можете упростить вычисления с "vecrminusvecs = ..." до "db = konstante ...", используя некоторую векторную алгебру; Вы делаете некоторые изменения масштаба, которые не кажутся необходимыми или, по крайней мере, могут быть оптимизированы для скорости.

edit : теперь я с подозрением отношусь к "норме"; sqrt(sum(x.^2,2)) работает с каждой строкой и, вероятно, быстрее, чем norm(), но вы должны измерить его, если хотите использовать самый быстрый подход.

Answer 2

Что вы хотите оптимизировать? Скорость? Первый совет по оптимизации:

Измерить.

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

Тем не менее, общее правило с Matlab - стараться избегать циклов (особенно вложенных) любой ценой и вместо этого выяснять, как представить ваши вычисления как матричные операции. Они быстрые и оптимизированные.

Answer 3

Я запустил этот код в профилировщике MATLAB, и несколько вещей выделяются:

1. Вы создаете и уничтожаете панель ожидания каждый раз вокруг цикла x = 1: xmax - вы должны держать панель ожидания открытой все время. Это заняло довольно много времени на моей машине.

2. Вам действительно нужно векторизовать хотя бы три внутренних цикла. Например, ваш расчет «dl» может быть заменен одним вызовом (что-то вроде - не проверено) dl = diff( leiter, 1, 1 ). Точно так же vecs = (leiter(1:N-1,:) + leiter(2:N,:))/2. Остальные выражения в этом цикле нужно немного доработать.

Answer 4

Это только ошибка на 20 стр / мин, но pi ~ = 3.1415927! = 3.1415297. (у вас есть опечатка, которая поменяла 2 и 9). Еще одна причина, кроме удобства использования встроенной константы.