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

2 次查看（过去 30 天）
SAHIL SAHOO
SAHIL SAHOO 2022-8-11
编辑： Torsten 2022-8-11
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

请先登录，再回答此问题。

采纳的回答

Alberto Cuadra Lara
编辑：Alberto Cuadra Lara 2022-8-11
Ran in:
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 个评论
显示 3更早的评论
SAHIL SAHOO
SAHIL SAHOO 2022-8-11
Ran in:
@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
Torsten 2022-8-11
编辑：Torsten 2022-8-11
Ran in:
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

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by