Plotting from a for loop

3 views (last 30 days)
JJH
JJH on 26 Jul 2017
Commented: dpb on 27 Jul 2017
Hi, I have a code (given below) where I calculate the trace distance between two matrices and take the output which I have called td2(t) where I have a for loop for values of t. I am able to make td2 into a matrix of all outputs of t and plot for when x and t go from 0 to 5. However when I try to change these values to a different amount, i.e. changing to x=0:1:100 (being sure to change this wherever it comes up in the code) I get an error saying the vectors must be the same length. I'm not sure which part isn't changing correctly. Any help would be appreciated!
A1=0.0023;
A6=12.5/3000;
we1=0;
we5=1.0674245;
r=0;
rhoE=[0.5 + 0.115047, 0.564418 - 0.0959956*1i; 0.564418 - 0.0959956*1i, 0.5 - 0.115047];
rhoN=[0 0; 0 1];
rhoEN=kron(rhoE,rhoN);
rhoSB=tp([0.5 0.5*r; 0.5*r 0.5],4);
rhoin=kron(rhoEN,rhoSB);
PauliX=[0 1; 1 0];
PauliY=[0 -1i; 1i 0];
PauliZ=[1 0; 0 -1];
td1=zeros(0:1:5);
for t=0:1:5;
H1=we1*kron(PauliX,eye(2))+A1*(kron(PauliX,PauliX)+kron(PauliY,PauliY)+kron(PauliZ,PauliZ));
H5=we5*megakron(PauliX,eye(2),eye(2),eye(2),eye(2),eye(2))+ ...
A6*(megakron(PauliX,PauliX,eye(2),eye(2),eye(2),eye(2))+ ...
megakron(PauliX,eye(2),PauliX,eye(2),eye(2),eye(2))+ ...
megakron(PauliX,eye(2),eye(2),PauliX,eye(2),eye(2))+ ...
megakron(PauliX,eye(2),eye(2),eye(2),PauliX,eye(2))+ ...
megakron(PauliX,eye(2),eye(2),eye(2),eye(2),PauliX)+ ...
megakron(PauliY,PauliY,eye(2),eye(2),eye(2),eye(2))+ ...
megakron(PauliY,eye(2),PauliY,eye(2),eye(2),eye(2))+ ...
megakron(PauliY,eye(2),eye(2),PauliY,eye(2),eye(2))+ ...
megakron(PauliY,eye(2),eye(2),eye(2),PauliY,eye(2))+ ...
megakron(PauliY,eye(2),eye(2),eye(2),eye(2),PauliY)+ ...
megakron(PauliZ,PauliZ,eye(2),eye(2),eye(2),eye(2))+ ...
megakron(PauliZ,eye(2),PauliZ,eye(2),eye(2),eye(2))+ ...
megakron(PauliZ,eye(2),eye(2),PauliZ,eye(2),eye(2))+ ...
megakron(PauliZ,eye(2),eye(2),eye(2),PauliZ,eye(2))+ ...
megakron(PauliZ,eye(2),eye(2),eye(2),eye(2),PauliZ));
U1=expm(-1i*H1*t);
U5=expm(-1i*H5*t);
rhoout=U5*rhoin*U5';
rhoENout=U1*rhoEN*U1';
pt=PartialTrace(rhoout,[3,4,5,6],[2,2,2,2,2,2]);
td1=(pt-rhoENout)^2;
td2(t+1)=sqrt(sum(td1(:)));
x=0:1:5;
plot(x,td2)
end
  2 Comments
Chad Greene
Chad Greene on 26 Jul 2017
Whoa, if you can simplify that code, you'll have a better chance of getting help. It's pretty tough to parse what's going on there. Try to remove anything you know does not contribute to the problem you're having.
dpb
dpb on 26 Jul 2017
I did it for you this time, but going forward select the code sections and then press the {}Code button to format code legibly...it's asinine that the forum uses word wrap as the default for a coding forum, granted, but that's the way it is...

Sign in to comment.

Accepted Answer

dpb
dpb on 26 Jul 2017
td1=zeros(0:1:5);
for t=0:1:5;
...
pt=PartialTrace(rhoout,[3,4,5,6],[2,2,2,2,2,2]);
td1=(pt-rhoENout)^2;
td2(t+1)=sqrt(sum(td1(:)));
x=0:1:5;
plot(x,td2)
end
There's no point in redefining x over and over as a constant...move it out of the loop. Then, having "magic numbers" buried in the code is error prone and painful to have to edit as you've already learned so make a parameter that only need to change it (or use a functional form instead of script and make it an argument of the function).
N=5; % upper limit on t vector
t=0:N; % define t
td1=zeros(size(t)); % preallocate output to match
for i=1:length(t) % iterate over elements of x
...
U1=expm(-1i*H1*t(i)); % now use element of t array...
U5=expm(-1i*H5*t(i));
...
pt=PartialTrace(rhoout,[3,4,5,6],[2,2,2,2,2,2]);
td1=(pt-rhoENout)^2;
td2(i)=sqrt(sum(td1(:))); % store output
end
plot(t,td2) % move plot until done then plot the arrays
The above uses 1:length of array as the indexing to avoid the "plus one" issue in avoiding the zero-th illegal array indexing and then plots the whole arrays when accumulated. This has the side effect of also allowing other than integer-valued t values if that's eventually desired.
On changing N to other sizes, the above should automagically take care of that, too...
  1 Comment
dpb
dpb on 27 Jul 2017
[JJH Answer moved to Comment -- dpb]
That's great, thank you!

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!