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

1 view (last 30 days)
SAHIL SAHOO
SAHIL SAHOO on 18 Dec 2022
Answered: Mathieu NOE on 19 Dec 2022
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

Sign in to answer this question.

Accepted Answer

Mathieu NOE
Mathieu NOE on 19 Dec 2022
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

More Answers (0)

Categories

Find more on Multirate Signal Processing 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!