Plotting from a for loop

2 次查看(过去 30 天)
JJH
JJH 2017-7-26
评论: dpb 2017-7-27
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 个评论
Chad Greene
Chad Greene 2017-7-26
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 2017-7-26
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...

请先登录,再进行评论。

采纳的回答

dpb
dpb 2017-7-26
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 个评论
dpb
dpb 2017-7-27
[JJH Answer moved to Comment -- dpb]
That's great, thank you!

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Annotations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by