how to add loop in this so that I can plot for 100 times and I will combine all in one plot.

1 次查看（过去 30 天）
SAHIL SAHOO
SAHIL SAHOO 2022-12-18
回答： Mathieu NOE 2022-12-19
ti = 0;
tf = 70E-8;
tspan=[ti tf];
KC = 1E-3;
y0= [(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      ((-3.14).*rand(20,1) + (3.14).*rand(20,1))];
yita_mn = [
    0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
    1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1;
    1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
    ]*(KC);
N = 20;
tp = 1E-12;
o = sort(10e2*rand(1,20),'ascend');
[T,Y]= ode45(@(t,y) rate_eq(t,y,yita_mn,N,o),tspan./tp,y0);
r = ((1/20).*( exp(i.*Y(:,3)) + exp(i.*Y(:,6)) + exp(i.*Y(:,9)) + exp(i.*Y(:,12)) + exp(i.*Y(:,15)) ...
    +exp(i.*Y(:,18)) +exp(i.*Y(:,21)) +exp(i.*Y(:,24)) + exp(i.*Y(:,27)) + exp(i.*Y(:,30)) + exp(i.*Y(:,33)) ...
    + exp(i.*Y(:,36)) + exp(i.*Y(:,39)) +exp(i.*Y(:,42)) + exp(i.*Y(:,45)) + exp(i.*Y(:,48)) + exp(i.*Y(:,51)) + exp(i.*Y(:,54))+ exp(i.*Y(:,57)) + exp(i.*Y(:,60))));
L = (((Y(:,2) - mean(Y(:,2))).^2 + (Y(:,5) - mean(Y(:,5))).^2 + (Y(:,8) - mean(Y(:,8))).^2 + (Y(:,11) - mean(Y(:,11))).^2 +(Y(:,14) - mean(Y(:,14))).^2 +(Y(:,17) - mean(Y(:,17))).^2+ (Y(:,20) - mean(Y(:,20))).^2+(Y(:,23) - mean(Y(:,23))).^2 + (Y(:,26) - mean(Y(:,26))).^2 +  (Y(:,29) - mean(Y(:,29))).^2+(Y(:,32) - mean(Y(:,32))).^2 + (Y(:,35) - mean(Y(:,35))).^2 + (Y(:,38) - mean(Y(:,38))).^2+(Y(:,41) - mean(Y(:,41))).^2 +(Y(:,44) - mean(Y(:,44))).^2 + (Y(:,47) - mean(Y(:,47))).^2+(Y(:,50) - mean(Y(:,50))).^2 +(Y(:,53) - mean(Y(:,53))).^2 + (Y(:,56) - mean(Y(:,56))).^2+(Y(:,59) - mean(Y(:,59))).^2)/19).^(0.5);
% I want to plot this for 100 times for different initial condition and combine all the plots in a single
% plot.
figure(1)
plot(T,L,'linewidth',1.5);
xlim([0 2E4]);
xlabel("time")
ylabel("standar deviation of amplitude")
grid on
set(gca,'fontname','times New Roman','fontsize',18,'linewidth',1.8);
figure(2)
plot(T,abs(r),'linewidth',1.5);
xlim([0 2E4]);
xlabel("time")
ylabel("orderparameter")
grid on
set(gca,'fontname','times New Roman','fontsize',18,'linewidth',1.8);
function dy = rate_eq(t,y,yita_mn,N,o)
dy = zeros(4*N,1);
dGdt = zeros(N,1);
dAdt = zeros(N,1);
dOdt = zeros(N,1);
P = 0.01;
a = 0.5;
T = 2000;
tp = 1E-12;
Gt = y(1:3:3*N-2);
At = y(2:3:3*N-1);
Ot = y(3:3:3*N-0);
k = 1E-3;
for i = 1:N
    dGdt(i) = (P - Gt(i) - (1 + 2.*Gt(i)).*((At(i)))^2)./T ;
    dAdt(i) = Gt(i)*(At(i));
    dOdt(i) = -a.*Gt(i) + o(1,i).*tp;  
    for j = 1:N
        dAdt(i) = dAdt(i) + yita_mn(i,j)*(At(j))*cos(Ot(j)-Ot(i));
        dOdt(i) = dOdt(i) + yita_mn(i,j)*((At(j)/At(i)))*sin(Ot(j)-Ot(i));
    end
end
dy(1:3:3*N-2) = dGdt;
dy(2:3:3*N-1) = dAdt;
dy(3:3:3*N-0) = dOdt;
n1 = (1:20)';
n2 = circshift(n1,-1);
n61 = n1 +60;
n62 = circshift(n61,-1);
n80 = circshift(n61,1);
j2 =  3*(1:20)-1;
j5 = circshift(j2,-1);
j8 = circshift(j2,-2);
j59 = circshift(j2,1);
dy(n61) = (o(1,n2).' - o(1,n1).').*tp - a.*(Gt(n2) - Gt(n1)) - (k).*(y(j2)./y(j5)).*sin(y(n61)) - (k).*(y( j5)./y(j2)).*sin(y(n61)) + (k).*(y(j8)./y(j5)).*sin(y(n62)) + (k).*(y(j59)./y(j2)).*sin(y(n80));  
end

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

采纳的回答

Mathieu NOE
Mathieu NOE 2022-12-19
hello
simply turn your main code into a function a call it 100 times in a for loop
clc
clearvars
close all
ti = 0;
tf = 70E-8;
tspan=[ti tf];
% I want to plot this for 100 times for different initial condition and combine all the plots in a single
% plot.
for ci = 1:100
[T,L,r] = main(tspan);
figure(1),hold on
plot(T,L,'linewidth',1.5);
xlim([0 2E4]);
xlabel("time")
ylabel("standar deviation of amplitude")
grid on
set(gca,'fontname','times New Roman','fontsize',18,'linewidth',1.8);
figure(2),hold on
plot(T,abs(r),'linewidth',1.5);
xlim([0 2E4]);
xlabel("time")
ylabel("orderparameter")
grid on
set(gca,'fontname','times New Roman','fontsize',18,'linewidth',1.8);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [T,L,r] = main(tspan)
KC = 1E-3;
y0= [(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(1e-2)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
((-3.14).*rand(20,1) + (3.14).*rand(20,1))];
yita_mn = [
0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1;
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
]*(KC);
N = 20;
tp = 1E-12;
o = sort(10e2*rand(1,N),'ascend');
[T,Y]= ode45(@(t,y) rate_eq(t,y,yita_mn,N,o),tspan./tp,y0);
r = ((1/20).*( exp(1i.*Y(:,3)) + exp(1i.*Y(:,6)) + exp(1i.*Y(:,9)) + exp(1i.*Y(:,12)) + exp(1i.*Y(:,15)) ...
+exp(1i.*Y(:,18)) +exp(1i.*Y(:,21)) +exp(1i.*Y(:,24)) + exp(1i.*Y(:,27)) + exp(1i.*Y(:,30)) + exp(1i.*Y(:,33)) ...
+ exp(1i.*Y(:,36)) + exp(1i.*Y(:,39)) +exp(1i.*Y(:,42)) + exp(1i.*Y(:,45)) + exp(1i.*Y(:,48)) + exp(1i.*Y(:,51)) + exp(1i.*Y(:,54))+ exp(1i.*Y(:,57)) + exp(1i.*Y(:,60))));
L = (((Y(:,2) - mean(Y(:,2))).^2 + (Y(:,5) - mean(Y(:,5))).^2 + (Y(:,8) - mean(Y(:,8))).^2 + (Y(:,11) - mean(Y(:,11))).^2 +(Y(:,14) - mean(Y(:,14))).^2 +(Y(:,17) - mean(Y(:,17))).^2+ (Y(:,20) - mean(Y(:,20))).^2+(Y(:,23) - mean(Y(:,23))).^2 + (Y(:,26) - mean(Y(:,26))).^2 + (Y(:,29) - mean(Y(:,29))).^2+(Y(:,32) - mean(Y(:,32))).^2 + (Y(:,35) - mean(Y(:,35))).^2 + (Y(:,38) - mean(Y(:,38))).^2+(Y(:,41) - mean(Y(:,41))).^2 +(Y(:,44) - mean(Y(:,44))).^2 + (Y(:,47) - mean(Y(:,47))).^2+(Y(:,50) - mean(Y(:,50))).^2 +(Y(:,53) - mean(Y(:,53))).^2 + (Y(:,56) - mean(Y(:,56))).^2+(Y(:,59) - mean(Y(:,59))).^2)/19).^(0.5);
end
function dy = rate_eq(~,y,yita_mn,N,o)
dy = zeros(4*N,1);
dGdt = zeros(N,1);
dAdt = zeros(N,1);
dOdt = zeros(N,1);
P = 0.01;
a = 0.5;
T = 2000;
tp = 1E-12;
Gt = y(1:3:3*N-2);
At = y(2:3:3*N-1);
Ot = y(3:3:3*N-0);
k = 1E-3;
for i = 1:N
dGdt(i) = (P - Gt(i) - (1 + 2.*Gt(i)).*((At(i)))^2)./T ;
dAdt(i) = Gt(i)*(At(i));
dOdt(i) = -a.*Gt(i) + o(1,i).*tp;
for j = 1:N
dAdt(i) = dAdt(i) + yita_mn(i,j)*(At(j))*cos(Ot(j)-Ot(i));
dOdt(i) = dOdt(i) + yita_mn(i,j)*((At(j)/At(i)))*sin(Ot(j)-Ot(i));
end
end
dy(1:3:3*N-2) = dGdt;
dy(2:3:3*N-1) = dAdt;
dy(3:3:3*N-0) = dOdt;
n1 = (1:20)';
n2 = circshift(n1,-1);
n61 = n1 +60;
n62 = circshift(n61,-1);
n80 = circshift(n61,1);
j2 = 3*(1:20)-1;
j5 = circshift(j2,-1);
j8 = circshift(j2,-2);
j59 = circshift(j2,1);
dy(n61) = (o(1,n2).' - o(1,n1).').*tp - a.*(Gt(n2) - Gt(n1)) - (k).*(y(j2)./y(j5)).*sin(y(n61)) - (k).*(y( j5)./y(j2)).*sin(y(n61)) + (k).*(y(j8)./y(j5)).*sin(y(n62)) + (k).*(y(j59)./y(j2)).*sin(y(n80));
end

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Signal Generation and Preprocessing 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by