my rand operation is applied correctly, I think but it's not working.
    4 views (last 30 days)
  
       Show older comments
    
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 Comments
Accepted Answer
  Chunru
      
      
 on 29 Aug 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 Comments
  Chunru
      
      
 on 29 Aug 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
More Answers (0)
See Also
Categories
				Find more on Develop Apps Using App Designer in Help Center and File Exchange
			
	Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!


