did I plot the "r" factor correctly in it?

1 view (last 30 days)
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?

Accepted Answer

Karim
Karim on 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

More Answers (0)

Categories

Find more on Historical Contests in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!