my rand operation is applied correctly, I think but it's not working.
3 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
SAHIL SAHOO
am 29 Aug. 2022
Kommentiert: Chunru
am 29 Aug. 2022
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 Kommentare
Akzeptierte Antwort
Chunru
am 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 Kommentare
Chunru
am 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
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Logical finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!