How to put matrices defined earlier in "for loop" later

I have this program defined,where E_i,E_j are 3x3 matrices, I use Qetlab Tensor Product to make 9x9 matrices

syms  x p_1 p_2 r
rho_alpha = (1/21).*[2 0 0 0 2 0 0 0 2;0 5-x 0 0 0 0 0 0 0;0 0 x 0 0 0 0 0 0;
    0 0 0 x 0 0 0 0 0;2 0 0 0 2 0 0 0 2;0 0 0 0 0 5-x 0 0 0;
    0 0 0 0 0 0 5-x 0 0;0 0 0 0 0 0 0 x 0;2 0 0 0 2 0 0 0 2];
E_1 = sqrt(r).*[1 0 0;0 sqrt(1-p_1) 0;0 0 sqrt(1-p_2)];
E_2 = sqrt(r).*[0 sqrt(p_1) 0;0 0 0;0 0 0];
E_3 = sqrt(r).*[0 0 sqrt(p_2);0 0 0;0 0 0];
E_4 = sqrt(1-r).*[sqrt(1-p_1-p_2) 0 0;0 1 0;0 0 1];
E_5 = sqrt(1-r).*[0 0 0; sqrt(p_1) 0 0;0 0 0];
E_6 = sqrt(1-r).*[0 0 0;0 0 0; sqrt(p_2) 0 0];
for i = 1:6
    for j = 1:6
        rho_alpha1 = Tensor(E_j,E_i)*rho_alpha*transpose(Tensor(E_i,E_j));
    end
end
rho_alpha1

the main issue i have is how to use E_1,E_2,E_3,E_4,E_5,E_6 matrices as E_i and E_j in for loop.

I am trying to replicate the results of this paper : Distillability sudden death and sudden birth in a

two-qutrit system under decoherence at finite temperature.

Answer 1

[Bruno Luong]

I have no idea if you expect this

syms  x p_1 p_2 r
rho_alpha = (1/21).*[2 0 0 0 2 0 0 0 2;0 5-x 0 0 0 0 0 0 0;0 0 x 0 0 0 0 0 0;
    0 0 0 x 0 0 0 0 0;2 0 0 0 2 0 0 0 2;0 0 0 0 0 5-x 0 0 0;
    0 0 0 0 0 0 5-x 0 0;0 0 0 0 0 0 0 x 0;2 0 0 0 2 0 0 0 2];
E = {sqrt(r).*[1 0 0;0 sqrt(1-p_1) 0;0 0 sqrt(1-p_2)];
     sqrt(r).*[0 sqrt(p_1) 0;0 0 0;0 0 0];
     sqrt(r).*[0 0 sqrt(p_2);0 0 0;0 0 0];
     sqrt(1-r).*[sqrt(1-p_1-p_2) 0 0;0 1 0;0 0 1];
     sqrt(1-r).*[0 0 0; sqrt(p_1) 0 0;0 0 0];
     sqrt(1-r).*[0 0 0;0 0 0; sqrt(p_2) 0 0]};
rho_alpha1 = 0;
n = length(E);
for i = 1:n
    for j = 1:n
        rho_alpha1 = rho_alpha1 + kron(E{j},E{i})*rho_alpha*transpose(kron(E{i},E{j}));
    end
end
rho_alpha1

rho_alpha1 = 

[Gurvir]

Hi bruno , thanks for helping out with that .

Could you tell me how you have written rho_alpha1 in this neat form ,where coefficents are defined.

I have almost completed the program, a little bit more help is required.

I have constructed a makeshift functional version of this program.

function rho_alpha1 = horodecki(x,p,r)
rho_alpha = (1/21).*[2 0 0 0 2 0 0 0 2;0 5-x 0 0 0 0 0 0 0;0 0 x 0 0 0 0 0 0;
    0 0 0 x 0 0 0 0 0;2 0 0 0 2 0 0 0 2;0 0 0 0 0 5-x 0 0 0;
    0 0 0 0 0 0 5-x 0 0;0 0 0 0 0 0 0 x 0;2 0 0 0 2 0 0 0 2];
E = {sqrt(r).*[1 0 0;0 sqrt(1-p) 0;0 0 sqrt(1-p)];
     sqrt(r).*[0 sqrt(p) 0;0 0 0;0 0 0];
     sqrt(r).*[0 0 sqrt(p);0 0 0;0 0 0];
     sqrt(1-r).*[sqrt(1-2*p) 0 0;0 1 0;0 0 1];
     sqrt(1-r).*[0 0 0; sqrt(p) 0 0;0 0 0];
     sqrt(1-r).*[0 0 0;0 0 0; sqrt(p) 0 0]};
rho_alpha1 = 0;
n = length(E);
for i = 1:n
    for j = 1:n
        rho_alpha1 = rho_alpha1 + kron(E{j},E{i})*rho_alpha*transpose(kron(E{i},E{j}));
    end
end
end

I use this program in constructing Two kinds of plots. The first program uses input from above function to generate data.

function data = negativity_data()
    arr = [];
    p_array = linspace(0, 0.5, 500);
    for p = p_array
        neg = Negativity(horodecki(4.3,p,0.9));
        arr = [arr; p neg];
    end
    data = arr;
end

Then i put it in another program to generate plot

function f = negativity_plot()
    d = negativity_data();
    x = d(:,1);
    y = d(:,2);
    f = plot(x,y);
end

Similarly ,i have another two programs ,

function data = realignment_data()
    arr = [];
    p_array = linspace(0, 0.5, 500);
    for p = p_array
        rea = 0.5*TraceNorm(Realignment(horodecki(4.3,p,0.9)))-0.5;
        arr = [arr; p rea];
    end
    data = arr;
end

%And the plot program

function f = realignment_plot()
    d = realignment_data();
    x = d(:,1);
    y = d(:,2);
    f = plot(x,y);
end

I am having two issues, (i) the plots are generated individually , i wanted both values in a single plot

(ii) whenever the values of the out hits 0, I want it to terminate,i.e I do not want plot where output values is negative.

[Bruno Luong]

Sorry I cannot help you. There seem to be no relationship between your original question and your plotting issue, and you post a bunch of functions (not programs) that are either unknown and not telling us how you call them, and I don't have symbolic toolbox (just run the code on Online server that has this toolbox).

Please open new question with MWE (Minimum Working Example). Please read the guideline how to ask question correctly in order to get fast answer.

[Gurvir]

Thanks for the reply. Will read the guideline and improve it.

But can you atleast tell me how you have written rho_alpha1 in a neat form ,where coefficents are defined.

[Bruno Luong]

I don't know what is "neat form"; Again I don't have symbolic toolbox and the knowledge to assist you more.

Please open a new thread if you have such similar question, and I hope you would clearly explain what mean "neat form"

[Gurvir]

I meant how do you get this from matlab

[Bruno Luong]

By this statement at the end of my code

rho_alpha1

the variable name (without semicolon), that displays the variable, don't you see it?

Again PLEASE POST A NEW THREAD.

Answer 2

[Voss]

Instead of making 6 separate variables, you can put those 6 matrices in a single variable that you can index, e.g., a 3D array or a cell array. Here they are in a 3D array:

syms  x p_1 p_2 r
rho_alpha = (1/21).*[2 0 0 0 2 0 0 0 2;0 5-x 0 0 0 0 0 0 0;0 0 x 0 0 0 0 0 0;
    0 0 0 x 0 0 0 0 0;2 0 0 0 2 0 0 0 2;0 0 0 0 0 5-x 0 0 0;
    0 0 0 0 0 0 5-x 0 0;0 0 0 0 0 0 0 x 0;2 0 0 0 2 0 0 0 2];
E = cat(3, ...
    sqrt(r).*[1 0 0;0 sqrt(1-p_1) 0;0 0 sqrt(1-p_2)], ...
    sqrt(r).*[0 sqrt(p_1) 0;0 0 0;0 0 0], ...
    sqrt(r).*[0 0 sqrt(p_2);0 0 0;0 0 0], ...
    sqrt(1-r).*[sqrt(1-p_1-p_2) 0 0;0 1 0;0 0 1], ...
    sqrt(1-r).*[0 0 0; sqrt(p_1) 0 0;0 0 0], ...
    sqrt(1-r).*[0 0 0;0 0 0; sqrt(p_2) 0 0] ...
    )

E(:,:,1) = 

E(:,:,2) = 

E(:,:,3) = 

E(:,:,4) = 

E(:,:,5) = 

E(:,:,6) = 

for i = 1:6
    for j = 1:6
        % note that rho_alpha1 is overwritten each time though the loops ...
        rho_alpha1 = Tensor(E(:,:,j),E(:,:,i))*rho_alpha*transpose(Tensor(E(:,:,i),E(:,:,j)));
    end
end

Unrecognized function or variable 'Tensor'.

rho_alpha1 % ... so this is rho_alpha1 for i == j == 6 only

[Gurvir]

thanks for the answer voss,but this is only giving rho_alpha1 for i==j==6 only.'kron' can be used instead of specialized package for 'Tensor'. But it is still giving only a single term. In a proper way there should be a sum over terms ,totalling 36 terms.

rho_alpha1 = kron(E(:,:,j),E(:,:,i))*rho_alpha*transpose(kron(E(:,:,i),E(:,:,j)));

Answer 3

[Walter Roberson]

syms  x p_1 p_2 r
rho_alpha = (1/21).*[2 0 0 0 2 0 0 0 2;0 5-x 0 0 0 0 0 0 0;0 0 x 0 0 0 0 0 0;
    0 0 0 x 0 0 0 0 0;2 0 0 0 2 0 0 0 2;0 0 0 0 0 5-x 0 0 0;
    0 0 0 0 0 0 5-x 0 0;0 0 0 0 0 0 0 x 0;2 0 0 0 2 0 0 0 2];
E{1} = sqrt(r).*[1 0 0;0 sqrt(1-p_1) 0;0 0 sqrt(1-p_2)];
E{2} = sqrt(r).*[0 sqrt(p_1) 0;0 0 0;0 0 0];
E{3} = sqrt(r).*[0 0 sqrt(p_2);0 0 0;0 0 0];
E{4} = sqrt(1-r).*[sqrt(1-p_1-p_2) 0 0;0 1 0;0 0 1];
E{5} = sqrt(1-r).*[0 0 0; sqrt(p_1) 0 0;0 0 0];
E{6} = sqrt(1-r).*[0 0 0;0 0 0; sqrt(p_2) 0 0];
for i = 1:6
    for j = 1:6
        rho_alpha = Tensor(E{j},E{i})*rho_alpha*transpose(Tensor(E{i},E{j}));
    end
end
rho_alpha

[Gurvir]

Hi walter,this program returns an error

E{1} = sqrt(r).*[1 0 0;0 sqrt(1-p_1) 0;0 0 sqrt(1-p_2)];

Invalid indexing assignment.

[Bruno Luong]

@Gurvir

clear E 

before running Walter code (always dangerous such assignment)

[Gurvir]

thanks bruno,yes the program now runs,but it has the same problem as that with voss's program. It only returns value for i == j == 6.It does not give a sum over the total 36 terms.

[Bruno Luong]

The flaw is not from Voss and Walter they simply duplicate your orginal code calculation inside the loops. Your original code doesn't do any sum, so how they suppose to know adding a sum on intermediate result?

[Gurvir]

Apologies,i am just in the learning phase of matlab.I thought the 'for' loops were for summing. Could you show me how to add sum in the program ?