Problem substituting a for loop with vectorization

I'm trying to speed up my code and I've been replacing a few for loops with arrays. However, I'm having trouble with this one and I'm not sure if it can be done the same way. The problem is that I need to calculate matrices, so when I introduce the arrays, instead of obtaining several matrices (one for each member of the array), I obtain only one giant matrix.

This is my current code:

for k1 = 0:(1*pi/180):(360*pi/180)
    for k2 = 0:(1*pi/180):(360*pi/180)
        
        JS = [
              -l1*cos(k1 + alfa1)        l1*sin(k1 + alfa1) ;
               l2*cos(k2 + alfa2)        l2*sin(k2 + alfa2)
             ];
     
        JE = [
              -sbr4.l1*cos(k1)           0                  ;
               0                        -sbr4.l2*cos(k2)
             ];

        js = det(JS);
        je = det(JE);
        
        if js == 0
            % (...)

        end

        if je == 0
            % (...)

        end
        
    end
end

And this is what I want:

k1 = linspace(0,360*pi/180,360);
k2 = linspace(0,360*pi/180,360);
K1 = repmat(k1,numel(k1),1);
K2 = repmat(transpose(k2),1,numel(k2));

JS = [
      -l1*cos(K1 + alfa1)        l1*sin(K1 + alfa1) ;
       l2*cos(K2 + alfa2)        l2*sin(K2 + alfa2)
     ];

JE = [
      -l1*cos(K1)                zeros(size(K1))    ;
       zeros(size(K1))          -l2*cos(K2)
     ];

js = det(JS);
je = det(JE);

if js == 0
    % (...)

end

if je == 0
    % (...)

end

The image to the left is what the first code gives me, which is what I want and the image to the right is what the second code gives me.

If someone shares how to solve it this way or how to speed it up using another method, I'd be grateful to know.

Answer 1

[Voss]

% just to have some values to run with:
l1 = 1;
l2 = 1;
alfa1 = 0;
alfa2 = 0;


k1 = linspace(0,2*pi,360);
k2 = linspace(0,2*pi,360);

% K1 = repmat(k1,numel(k2),1);
% K2 = repmat(transpose(k2),1,numel(k1));
[K1,K2] = meshgrid(k1,k2);

% K1 and K2 are 360-by-360.
% Make them 1-by-1-by-360^2:
K1 = shiftdim(K1(:),-2);
K2 = shiftdim(K2(:),-2);

% then JS and JE are 2-by-2-by-360^2

JS = [
      -l1*cos(K1 + alfa1)        l1*sin(K1 + alfa1) ;
       l2*cos(K2 + alfa2)        l2*sin(K2 + alfa2)
     ];

JE = [
      -l1*cos(K1)                zeros(size(K1))    ;
       zeros(size(K1))          -l2*cos(K2)
     ];

% you can calculate the determinant of each 2-by-2 matrix in JS and JE in
% a loop:
for ii = 1:size(JS,3)
    
    js = det(JS(:,:,ii));
    je = det(JE(:,:,ii));
    
    if js == 0
        % (...)
    end
    
    if je == 0
        % (...)
    end
    
end

[Kerman Bilbao Pinedo]

Thanks for your answer, but I have a new problem now.

I didn't mention it on the original post because I thought it wouldn't matter, but in the lines I wrote "(...)" inside the if statements I plot some variables that I obtain via K1 and K2 and have the same size as them. However, now I get an error because I can't plot more than 2 dimensions. I'm aware of the existence of 3D plots, but as I don't really understand what you are doing with the sizes of K1 and K2, I can't figure out how to solve this.

Here are the lines that give me the error:

if js == 0
    [Xp,Yp] = Delta_sbr4_resolucion_fi_gdl(sbr4,K1,K2); % This is just the script I use to obtain Xp and Yp via K1 and K2.
    plot(Xp,Yp,'.b')
    
end
    
if je == 0
    [Xp,Yp] = Delta_sbr4_resolucion_fi_gdl(sbr4,K1,K2);
    plot(Xp,Yp,'.r')
    
end

And this is the error message:

Error using plot
Data cannot have more than 2 dimensions.

Error in Delta_sbr4_main (line 171)
            plot(Xp,Yp,'.r')

Also, I should've mentioned before that I only want to plot those values that come from the values in K1 and K2 that make the determinants zero.

[Torsten]

At least you should tell us the sizes of Xp and Yp ...

[Kerman Bilbao Pinedo]

They have the same sizes as K1 and K2.

[Torsten]

Then

plot(squeeze(Xp),squeeze(Yp))

should work.

[Voss]

@Kerman Bilbao Pinedo: Try this:

% you can calculate the determinant of each 2-by-2 matrix in JS and JE in
% a loop:
for ii = 1:size(JS,3)
    
    js = det(JS(:,:,ii));
    je = det(JE(:,:,ii));
    
    if js == 0
        [Xp,Yp] = Delta_sbr4_resolucion_fi_gdl(sbr4,K1(ii),K2(ii)); % This is just the script I use to obtain Xp and Yp via K1 and K2.
        plot(Xp,Yp,'.b')
    end
    
    if je == 0
        [Xp,Yp] = Delta_sbr4_resolucion_fi_gdl(sbr4,K1(ii),K2(ii));
        plot(Xp,Yp,'.r')
    end
    
end

[Kerman Bilbao Pinedo]

Thanks a lot, it works!

[Voss]

You're welcome!