why the plot is not coming, is my plot function is wrong? please find the defect in the code.

1 view (last 30 days)
ti = 0;
tf = 1E-3;
tspan=[ti tf];
y0=[1; 1; 0; 1; 1; 0; 1; 1; 0; 1; 1; 0; 1; 1; 0]*10E-2;
yita_mn = [0 1 0 0 1; 1 0 1 0 0; 0 1 0 1 0; 0 0 1 0 1; 1 0 0 1 0]*(0.01);
N = 5;
tp = 5.4E-9;
[T,Y]= ode45(@(t,y) rate_eq(t,y,yita_mn,N),tspan,y0);
subplot 211
plot(T./tp,Y(:,2));
xlabel('time')
ylabel('amplitude')
subplot 212
plot(T./tp, mean(exp(Y(:,3)) + exp(Y(:,6))+ exp(Y(:,9)) +exp(Y(:,12))+exp(Y(:,15))));
xlabel('time')
ylabel('phase')
function dy = rate_eq(t,y,yita_mn,N)
dy = zeros(3*N,1);
dNdt = zeros(N,1);
dSdt = zeros(N,1);
dWdt = zeros(N,1);
P = 16.7;
a = 0.1;
tc = 230E-6;
tp =5.4E-9;
o = 3000;
k_c = 1/tp;
Nt = y(1:3:3*N-2);
St = y(2:3:3*N-1);
Wt = y(3:3:3*N-0);
for i = 1:N
dNdt(i) = (P - Nt(i).*((abs(St(i)))^2 +1))./tc ;
dSdt(i) = (Nt(i)-a).*((St(i))./tp);
dWdt(i) = o;
for j = 1:N
dSdt(i) = dSdt(i)+yita_mn(i,j)*k_c*(St(j))*cos(Wt(j)-Wt(i));
dWdt(i) = dWdt(i)+yita_mn(i,j)*k_c*((St(j)/St(i)))*sin(Wt(j)-Wt(i));
end
end
dy(1:3:3*N-2) = dNdt;
dy(2:3:3*N-1) = dSdt;
dy(3:3:3*N-0) = dWdt;
end

Accepted Answer

Alberto Cuadra Lara
Alberto Cuadra Lara on 11 Aug 2022
Edited: Alberto Cuadra Lara on 11 Aug 2022
Hello Sahil,
You are plotting a point not a vector with the same length as the time variable. Remove the mean function. I also suggest also using the log scale for time.
Here you have:
ti = 0;
tf = 1E-3;
tspan=[ti tf];
y0=[1; 1; 0; 1; 1; 0; 1; 1; 0; 1; 1; 0; 1; 1; 0]*10E-2;
yita_mn = [0 1 0 0 1; 1 0 1 0 0; 0 1 0 1 0; 0 0 1 0 1; 1 0 0 1 0]*(0.01);
N = 5;
tp = 5.4E-9;
[T,Y]= ode45(@(t,y) rate_eq(t,y,yita_mn,N),tspan,y0);
subplot 211
plot(T./tp, Y(:,2));
xlabel('time')
ylabel('amplitude')
set(gca, 'XScale', 'log')
subplot 212
plot(T./tp, exp(Y(:,3)) + exp(Y(:,6)) + exp(Y(:,9)) + exp(Y(:,12)) + exp(Y(:,15)));
xlabel('time')
ylabel('phase')
set(gca, 'XScale', 'log')
function dy = rate_eq(t, y, yita_mn, N)
dy = zeros(3*N,1);
P = 16.7;
a = 0.1;
tc = 230E-6;
tp =5.4E-9;
o = 3000;
k_c = 1/tp;
Nt = y(1:3:3*N-2);
St = y(2:3:3*N-1);
Wt = y(3:3:3*N-0);
for i = N:-1:1
dNdt(i) = (P - Nt(i) .* (St(i).^2 + 1)) ./ tc ;
dSdt(i) = (Nt(i) - a) .* St(i) ./ tp;
dWdt(i) = o;
for j = N:-1:1
dSdt(i) = dSdt(i) + yita_mn(i,j) .* k_c .* St(j) .* cos(Wt(j) - Wt(i));
dWdt(i) = dWdt(i) + yita_mn(i,j) .* k_c .* St(j) ./ St(i) .* sin(Wt(j) - Wt(i));
end
end
dy(1:3:3*N-2) = dNdt;
dy(2:3:3*N-1) = dSdt;
dy(3:3:3*N-0) = dWdt;
end
  5 Comments
Torsten
Torsten on 11 Aug 2022
Edited: Torsten on 11 Aug 2022
In the meantime, there were so many cases where you assigned a vector to a scalar. The same here. Don't you learn from the corrections the forum supplies ?
tspan= 0: 0.1E-3 :1e-3;
y0=[1; 1; 0; 1; 1; 0; 1; 1; 0; 1; 1; 0; 1; 1; 0]*10E-2;
yita_mn = [0 1 0 0 1; 1 0 1 0 0; 0 1 0 1 0; 0 0 1 0 1; 1 0 0 1 0]*(0.001);
N = 5;
tp = 5.4E-9;
[T,Y]= ode45(@(t,y) rate_eq(t,y,yita_mn,N),tspan,y0);
%A = zeros(length(tspan),1) ;
%for i = 1:length(tspan)
A = ((((Y(:,2) - mean(Y(:,2))).^2)).^(0.5))./(Y(:,2));
%end
r = zeros(length(tspan),1) ;
%for i = 1:length(tspan)
r = (exp(Y(:,3)) + exp(Y(:,6)) + exp(Y(:,9)) + exp(Y(:,12)) + exp(Y(:,15))) ;
%end
subplot 211
plot(T./tp,A);
xlabel('time')
ylabel('amplitude spread')
subplot 212
plot(T./tp,r);
xlabel('time')
ylabel('order parameter')
function dy = rate_eq(t,y,yita_mn,N)
dy = zeros(3*N,1);
dNdt = zeros(N,1);
dSdt = zeros(N,1);
dWdt = zeros(N,1);
P = 16.7;
a = 0.1;
tc = 230E-6;
tp =5.4E-9;
o = 20000;
k_c = 1/tp;
Nt = y(1:3:3*N-2);
St = y(2:3:3*N-1);
Wt = y(3:3:3*N-0);
for i = 1:N
dNdt(i) = (P - Nt(i).*((abs(St(i)))^2 +1))./tc ;
dSdt(i) = (Nt(i)-a).*((St(i))./tp);
dWdt(i) = o;
for j = 1:N
dSdt(i) = dSdt(i)+yita_mn(i,j)*k_c*(St(j))*cos(Wt(j)-Wt(i));
dWdt(i) = dWdt(i)+yita_mn(i,j)*k_c*((St(j)/St(i)))*sin(Wt(j)-Wt(i));
end
end
dy(1:3:3*N-2) = dNdt;
dy(2:3:3*N-1) = dSdt;
dy(3:3:3*N-0) = dWdt;
end

Sign in to comment.

More Answers (0)

Categories

Find more on Video games in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!