my rand operation is applied correctly, I think but it's not working.
    4 vues (au cours des 30 derniers jours)
  
       Afficher commentaires plus anciens
    
clear
ti = 0;
tf = 10E-3;
tspan=[ti tf];
y0=[1; 1; 0; 1; 1; 0; 1; 1; 0; 1; 1; 0; 1; 1; 0; 0; 0 ; 0; 0; 0]*10E-4;
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]*(1E-3);
N = 5;
tp = 5.4E-9;
[T,Y]= ode45(@(t,y) rate_eq(t,y,yita_mn,N),tspan,y0);
% time =T./tp done, because to measure the time in unit of tp 
plot(T./tp,Y(:,16));
function dy = rate_eq(t,y,yita_mn,N)
dy = zeros(4*N,1);
dGdt = zeros(N,1);
dAdt = zeros(N,1);
dOdt = zeros(N,1);
P = 16.7;
a = 0.1;
tc = 230E-6;
tp =5.4E-9;
o = [1 2 32 421 2]*10E2;
k_c = 1/tp;
Gt = y(1:3:3*N-2);
At = y(2:3:3*N-1);
Ot = y(3:3:3*N-0);
%A is the amplitude
%G is the gain
%O is the phase
for i = 1:N
    dGdt(i) = (P - Gt(i).*((abs(At(i)))^2 +1))./tc ; 
    dAdt(i) = (Gt(i)-a).*((At(i))./tp);
    dOdt(i) = o(1,:);  
    for j = 1:N
        dAdt(i) = dAdt(i)+yita_mn(i,j)*k_c*(At(j))*cos(Ot(j)-Ot(i));
        dOdt(i) = dOdt(i)+yita_mn(i,j)*k_c*((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;
dy(16)= dy(6) - dy(3);
dy(17)= dy(9) - dy(6);
dy(18)= dy(12) - dy(6);
dy(19)= dy(15) - dy(12);
dy(20)= dy(3) - dy(15);
end
0 commentaires
Réponse acceptée
  Chunru
      
      
 le 29 Août 2022
        Check out the line: 
dOdt(i) = o(i);      %o(1,:); 
ti = 0;
tf = 10E-3;
tspan=[ti tf];
y0=[1; 1; 0; 1; 1; 0; 1; 1; 0; 1; 1; 0; 1; 1; 0; 0; 0 ; 0; 0; 0]*10E-4;
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]*(1E-3);
N = 5;
tp = 5.4E-9;
[T,Y]= ode45(@(t,y) rate_eq(t,y,yita_mn,N),tspan,y0);
plot(T./tp,Y(:,16));
function dy = rate_eq(t,y,yita_mn,N)
dy = zeros(4*N,1);
dGdt = zeros(N,1);
dAdt = zeros(N,1);
dOdt = zeros(N,1);
P = 16.7;
a = 0.1;
tc = 230E-6;
tp =5.4E-9;
o = [1 2 32 421 2]*10E2;
k_c = 1/tp;
Gt = y(1:3:3*N-2);
At = y(2:3:3*N-1);
Ot = y(3:3:3*N-0);
%A is the amplitude
%G is the gain
%O is the phase
for i = 1:N
    dGdt(i) = (P - Gt(i).*((abs(At(i)))^2 +1))./tc ; 
    dAdt(i) = (Gt(i)-a).*((At(i))./tp);
    dOdt(i) = o(i);      %o(1,:);  
    for j = 1:N
        dAdt(i) = dAdt(i)+yita_mn(i,j)*k_c*(At(j))*cos(Ot(j)-Ot(i));
        dOdt(i) = dOdt(i)+yita_mn(i,j)*k_c*((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;
dy(16)= dy(6) - dy(3);
dy(17)= dy(9) - dy(6);
dy(18)= dy(12) - dy(6);
dy(19)= dy(15) - dy(12);
dy(20)= dy(3) - dy(15);
end
2 commentaires
  Chunru
      
      
 le 29 Août 2022
				It seems that "o = rand(1,5)*10E2" should be a variable outside of function rate_eq.
ti = 0;
tf = 10E-3;
tspan=[ti tf];
y0=[1; 1; 0; 1; 1; 0; 1; 1; 0; 1; 1; 0; 1; 1; 0; 0; 0 ; 0; 0; 0]*10E-4;
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]*(1E-3);
N = 5;
o = rand(1,5)*10E2   % <==============
tp = 5.4E-9;
[T,Y]= ode45(@(t,y) rate_eq(t,y,yita_mn,N, o),tspan,y0);
plot(T./tp,Y(:,20));
function dy = rate_eq(t,y,yita_mn,N, o)   % o as parameters
dy = zeros(4*N,1);
dGdt = zeros(N,1);
dAdt = zeros(N,1);
dOdt = zeros(N,1);
P = 16.7;
a = 0.1;
tc = 230E-6;
tp =5.4E-9;
k_c = 1/tp;
Gt = y(1:3:3*N-2);
At = y(2:3:3*N-1);
Ot = y(3:3:3*N-0);
%A is the amplitude
%G is the gain
%O is the phase
for i = 1:N
    dGdt(i) = (P - Gt(i).*((abs(At(i)))^2 +1))./tc ; 
    dAdt(i) = (Gt(i)-a).*((At(i))./tp);
    dOdt(i) = o(1,i);      %o(1,:);  
    for j = 1:N
        dAdt(i) = dAdt(i)+yita_mn(i,j)*k_c*(At(j))*cos(Ot(j)-Ot(i));
        dOdt(i) = dOdt(i)+yita_mn(i,j)*k_c*((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;
dy(16)= dy(6) - dy(3);
dy(17)= dy(9) - dy(6);
dy(18)= dy(12) - dy(6);
dy(19)= dy(15) - dy(12);
dy(20)= dy(3) - dy(15);
end
Plus de réponses (0)
Voir également
Catégories
				En savoir plus sur Develop Apps Using App Designer dans Help Center et File Exchange
			
	Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!



