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

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

Answer 1

[Alberto Cuadra Lara]

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

[SAHIL SAHOO]

but I want to plot the mean of (exp(Y(:,3)) + exp(Y(:,6)) + exp(Y(:,9)) + exp(Y(:,12)) + exp(Y(:,15))) for the different time. because I get Y(:,3) Y(:,6) Y(:,9) Y(:,12) and Y(:,15) for different time, how can i do that?

[SAHIL SAHOO]

@Alberto Cuadra Lara is this correct

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(i) = ((mean((Y(:,2) - mean(Y(:,2))).^2)).^(0.5))./mean(Y(:,2));
end

r = zeros(length(tspan),1) ; 
for i = 1:length(tspan)
    r(i) = mean(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

[Walter Roberson]

Constant output does not seem the most desirable to me ?

[SAHIL SAHOO]

@Walter Roberson I think this code should be right, still I get error. could you please check what's wrong in it.

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(i) = ((((Y(:,2) - mean(Y(:,2))).^2)).^(0.5))./(Y(:,2));
end

Unable to perform assignment because the left and right sides have a different number of elements.

r = zeros(length(tspan),1) ; 
for i = 1:length(tspan)
     r(i) = (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

[Torsten]

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