Double for loop to generate vector from sequential results

Hello,

I am trying to generate a "Double for loop", the general idea is execute some operations with matrices, but there is a scalar variable multiplying some of the matrices. Hence, I want the variable "om" to change from 1.01 to 1.99 at 0.01 increments, next calculate the eigenvalues of the resulting matrix, and finally store them in a vector array from 1 to 99. But, what happens in my code is that matlab solves all "om" values, and then stores the very last calculated eigenvalue 99 times.

This is the code thus far:

%Define Matrices
A1 = [52 -1 0 -51;...
      -51 52 -1 0;...
      0 -51 52 -1;...
      -1 0 -51 52];
 
A2 = [52 -1 0 -0;...
      -51 52 -1 0;...
      0 -51 52 -1;...
      0 0 -51 52];
 
%L-U Decomposition
[L1,U1,P1] = lu(A1); %Matrix A1
D1 = diag(diag(A1));
[L2,U2,P2] = lu(A2); %Matrix A2
D2 = diag(diag(A2));
 
%Iteration matrix (M1) for A1
 
vs = 1.01:0.01:1.99;
s = size(vs,2);
ev1 = 1:s;
eM1 = 1:s;
 
for i = 1:s
    for om = 1.01:0.01:1.99
     M1 = inv(D1-om*L1)*((1-om)*D1+om*U1); %Generate M1 varying omega
     eM1 = eig(M1); %Calculate eigenvalues
     %e1M1 = eM1(1,1);
     %e2M1 = eM1(2,1);
     %e3M1 = eM1(3,1);
     %e4M1 = eM1(4,1);
     ev1(i) = eM1(1,1);
    end
end

The A2, and D2 values currently are not used.

Thank you

Answer 1

[Stephen23]

You are overthinking things, you only need one loop:

vs  = 1.01:0.01:1.99;
s   = size(vs,2);
ev1 = nan(1,s); % preallocate
for k = 1:s
    om = vs(k);
    M1 = inv(D1-om*L1)*((1-om)*D1+om*U1); %Generate M1 varying omega
    eM1 = eig(M1); %Calculate eigenvalues
    %e1M1 = eM1(1,1);
    %e2M1 = eM1(2,1);
    %e3M1 = eM1(3,1);
    %e4M1 = eM1(4,1);
    ev1(k) = eM1(1,1);
end

Note that inv is not recommended for solving systems of equations: read the inv help for info on how to replace it with more efficient and numerically precise operations. You should use mldivide:

(D1-om*L1)\((1-om)*D1+om*U1)
%         ^ Read the INV documentation!

[Manuel Umanzor]

Thanks a lot

[Christine Tobler]

It looks like you're effectively computing eig(inv(B)*A), with A = ((1-om)*D1+om*U1) and B = D1-om*L1. A more accurate method is to compute eig(A, B) instead, which solves the system A*x = lambda*B*x (as opposed to the current method of solving inv(B)*A*x = lambda*x).