did I plot the "r" factor correctly in it?
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
SAHIL SAHOO
am 27 Aug. 2022
Beantwortet: Karim
am 29 Aug. 2022
ti = 0;
tf = 1E-4;
tspan=[ti tf];
y0=[1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 0; 0; 0; 0; 0];
N = 5;
tp = 5.4E-9;
[T,Y]= ode45(@(t,y) rate_eq(t,y,N),tspan,y0);
% did this r part is correct?
r = (exp((0+1i).*Y(:,11)) + exp((0+1i).*Y(:,12)) + exp((0+1i).*Y(:,13)) + exp((0+1i).*Y(:,14)) + exp((0+1i).*Y(:,15)))./5 ;
plot(T./tp,r);
xlabel('time')
ylabel('order parameter')
function dy = rate_eq(t,y,N)
dy = zeros(3*N,1);
P = 16.7;
a = 0.1;
tc = 230E-6;
tp = 5.4E-9;
o1 = 5308;
o2 = 4008;
o3 = 7009;
o4 = 6008;
o5 = 7080;
k1 = 1E-3;
k2 = 1E-3;
k3 = 1E-3;
k4 = 1E-3;
k5 = 1E-3;
dy(1) = (P - y(1).*((abs(y(6)))^2 +1))./tc;
dy(2) = (P - y(2).*((abs(y(7)))^2 +1))./tc;
dy(3) = (P - y(3).*((abs(y(8)))^2 +1))./tc;
dy(4) = (P - y(4).*((abs(y(9)))^2 +1))./tc;
dy(5) = (P - y(5).*((abs(y(10)))^2 +1))./tc;
dy(6)= (y(1)-a).*((y(6))./tp) + (k1./tp).*(y(7)).*cos(y(11)) + (k5./tp).*(y(10))*cos(y(15));
dy(7)= (y(2)-a).*((y(7))./tp) + (k2./tp).*(y(8)).*cos(y(12)) + (k1./tp).*(y(6))*cos(y(11));
dy(8)= (y(3)-a).*((y(8))./tp) + (k3./tp).*(y(9)).*cos(y(13)) + (k2./tp).*(y(7))*cos(y(12));
dy(9)= (y(4)-a).*((y(9))./tp) + (k4./tp).*(y(10)).*cos(y(14)) + (k3./tp).*(y(8))*cos(y(13));
dy(10)= (y(5)-a).*((y(10))./tp) + (k5./tp).*(y(9)).*cos(y(14)) +(k4./tp).* (y(6))*cos(y(15));
dy(11) = o1 - (k1./tp).*(y(6)/y(7)).*(sin(y(11))) - (k1./tp).*(y(10)/y(6)).*(sin(y(11))) + (k2./tp).*(y(8)./y(7))*sin(y(12)) + (k5./tp).*(y(10)/y(6)).*sin(y(15));
dy(12) = o2 - (k2./tp).*(y(7)/y(8)).*(sin(y(12))) - (k2./tp).*(y(8)/y(7)).*(sin(y(12))) + (k3./tp).*(y(9)./y(8))*sin(y(13)) + (k1./tp).*(y(6)/y(7)).*sin(y(11));
dy(13) = o3 - (k3./tp).*(y(8)/y(9)).*(sin(y(13))) - (k3./tp).*(y(9)/y(8)).*(sin(y(13))) + (k4./tp).*(y(10)./y(9))*sin(y(14)) + (k2./tp).*(y(7)/y(8)).*sin(y(12));
dy(14) = o4 - (k4./tp).*(y(9)/y(10)).*(sin(y(14))) - (k4./tp).*(y(10)/y(9)).*(sin(y(14))) + (k5./tp).*(y(6)./y(10))*sin(y(15)) + (k3./tp).*(y(8)/y(9)).*sin(y(13));
dy(15) = o5 - (k5./tp).*(y(10)/y(6)).*(sin(y(15))) - (k5./tp).*(y(6)/y(10)).*(sin(y(15))) + (k1./tp).*(y(7)./y(6))*sin(y(11)) + (k4./tp).*(y(9)/y(10)).*sin(y(14));
end
============================
= r, this is the equation of r here and i want to calculate the r value, and the ϕ1=Y(:,11),ϕ2 = Y(:,12),ϕ3 = Y(:,13), ϕ4 = Y(:,14), ϕ5 = Y(:,15) and N =5
please tell me did I plot r correctly or not?
0 Kommentare
Akzeptierte Antwort
Karim
am 29 Aug. 2022
at first glance you still need to take the absolute value, see below
ti = 0;
tf = 1E-4;
tspan=[ti tf];
y0=[1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 0; 0; 0; 0; 0];
N = 5;
tp = 5.4E-9;
[T,Y]= ode45(@(t,y) rate_eq(t,y,N),tspan,y0);
% evaluate r
r = ( exp(1i*Y(:,11)) + exp(1i*Y(:,12)) + exp(1i*Y(:,13)) + exp(1i*Y(:,14)) + exp(1i*Y(:,15)) ) / N;
r = abs( r );
plot(T./tp,r);
xlabel('time')
ylabel('order parameter')
function dy = rate_eq(t,y,N)
dy = zeros(3*N,1);
P = 16.7;
a = 0.1;
tc = 230E-6;
tp = 5.4E-9;
o1 = 5308;
o2 = 4008;
o3 = 7009;
o4 = 6008;
o5 = 7080;
k1 = 1E-3;
k2 = 1E-3;
k3 = 1E-3;
k4 = 1E-3;
k5 = 1E-3;
dy(1) = (P - y(1).*((abs(y(6)))^2 +1))./tc;
dy(2) = (P - y(2).*((abs(y(7)))^2 +1))./tc;
dy(3) = (P - y(3).*((abs(y(8)))^2 +1))./tc;
dy(4) = (P - y(4).*((abs(y(9)))^2 +1))./tc;
dy(5) = (P - y(5).*((abs(y(10)))^2 +1))./tc;
dy(6)= (y(1)-a).*((y(6))./tp) + (k1./tp).*(y(7)).*cos(y(11)) + (k5./tp).*(y(10))*cos(y(15));
dy(7)= (y(2)-a).*((y(7))./tp) + (k2./tp).*(y(8)).*cos(y(12)) + (k1./tp).*(y(6))*cos(y(11));
dy(8)= (y(3)-a).*((y(8))./tp) + (k3./tp).*(y(9)).*cos(y(13)) + (k2./tp).*(y(7))*cos(y(12));
dy(9)= (y(4)-a).*((y(9))./tp) + (k4./tp).*(y(10)).*cos(y(14)) + (k3./tp).*(y(8))*cos(y(13));
dy(10)= (y(5)-a).*((y(10))./tp) + (k5./tp).*(y(9)).*cos(y(14)) +(k4./tp).* (y(6))*cos(y(15));
dy(11) = o1 - (k1./tp).*(y(6)/y(7)).*(sin(y(11))) - (k1./tp).*(y(10)/y(6)).*(sin(y(11))) + (k2./tp).*(y(8)./y(7))*sin(y(12)) + (k5./tp).*(y(10)/y(6)).*sin(y(15));
dy(12) = o2 - (k2./tp).*(y(7)/y(8)).*(sin(y(12))) - (k2./tp).*(y(8)/y(7)).*(sin(y(12))) + (k3./tp).*(y(9)./y(8))*sin(y(13)) + (k1./tp).*(y(6)/y(7)).*sin(y(11));
dy(13) = o3 - (k3./tp).*(y(8)/y(9)).*(sin(y(13))) - (k3./tp).*(y(9)/y(8)).*(sin(y(13))) + (k4./tp).*(y(10)./y(9))*sin(y(14)) + (k2./tp).*(y(7)/y(8)).*sin(y(12));
dy(14) = o4 - (k4./tp).*(y(9)/y(10)).*(sin(y(14))) - (k4./tp).*(y(10)/y(9)).*(sin(y(14))) + (k5./tp).*(y(6)./y(10))*sin(y(15)) + (k3./tp).*(y(8)/y(9)).*sin(y(13));
dy(15) = o5 - (k5./tp).*(y(10)/y(6)).*(sin(y(15))) - (k5./tp).*(y(6)/y(10)).*(sin(y(15))) + (k1./tp).*(y(7)./y(6))*sin(y(11)) + (k4./tp).*(y(9)/y(10)).*sin(y(14));
end
0 Kommentare
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Historical Contests finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!