Matlab: Magnetic fields calculated with biot-savarts-law

Question

Hello! This is the code for the computation of magnetic fields using the biot-savart-law. I hope to get some tipps for the optimization of this code. Regrettably I use german language :( I never will do this again. :)

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.1415297;
spulenlaenge = 20; % Spulenlaenge in mm
steigung = spulenlaenge / windungen
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
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

What do you want to optimize? Speed? The first tip of optimization:

Measure it.

So, first measure where you spend most of the time (probably in some loop). Put timers around so that you get a good feel of it. Then see if you could do something about it. You can't optimize as long as you don't know what you need to optimize.

That said, the general rule with Matlab is to try to avoid loops (especially nested ones) at all costs and instead figure out how to present your computations as matrix operations. They are fast and optimized.

Answer 2

not necessarily speed optimization, but some notes:

• use pi, not 3.14159 or 3.1415927
• you can generally vectorize your loops:

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

vectorized:

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

Most matlab functions will take vectors/matrices as arguments, including cos(), sin(), exp(), log(), etc.

For low numbers of elements (say < a few hundred) it may not be worth the effort to vectorize.

Vector magnitude: instead of sqrt(dB(1).^2 + dB(2).^2 + dB(3).^2) use norm(dB) (note that norm does NOT operate on a matrix in a row-wise fashion but rather on the whole) though that won't save much

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

consider changing to

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

Why are you calculating r using a square root when you are just squaring it later?

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

Change to

r2 = sum(vecrminusvecs.^2);

and use r2 in place of r.^2

My guess is you can probably simplify the calculations from "vecrminusvecs = ..." to "db = konstante..." by using some vector algebra; you're doing some rescaling that doesn't quite seem like it's necessary, or at least could be optimized for speed.

edit: I am now suspicious of "norm"; sqrt(sum(x.^2,2)) operates on each row and is probably faster than norm() but you should measure it if you want to use the fastest approach.

Answer 3

I ran this code in the MATLAB profiler, and a couple of things stand out:

1. You are creating and destroying the waitbar each time around the x = 1:xmax loop - you should keep the waitbar open the whole time. This was taking a fair amount of time on my machine.

2. You really do need to vectorize at least the three inner loops. For example, your calculation of "dl" could be replaced by a single call to (something like - untested) dl = diff( leiter, 1, 1 ). Likewise, vecs = (leiter(1:N-1,:) + leiter(2:N,:))/2. The remaining expressions in that loop need a bit more work to tease apart.

Answer 4

This is only a 20ppm error, but pi ~= 3.1415927 != 3.1415297. (you have a typo which switched the 2 and the 9). Another reason besides convenience to use the built-in constant.