in this I apply a loop so that I will have different value of Y(:,11) for the different value of initial condition for Y(:,11), what is correct algorithm for this??
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
SAHIL SAHOO
am 22 Aug. 2022
Beantwortet: Sam Chak
am 22 Aug. 2022
please check this code is it right or worng? and please lemme know he correct program for this.
ti = 0;
tf = 1E-3;
tspan=[ti tf];
k = 1:5000:2*pi;
for i = 1:length(k)
y0(i)=[1; 1; 1; 1; 1; 1; 1; 1; 1; 1; k(i); 0; 0; 0; 0];
N = 5;
tp = 5.4E-9;
[T,Y]= ode45(@(t,y) rate_eq(t,y,N),tspan,y0(i));
E(i) = mean(Y(:,11));
q(i) = (5.*(log(E(i))))./(2*pi);
end
q
function dy = rate_eq(t,y,N)
dy = zeros(3*N,1);
P = 1.67;
a = 0.1;
tc = 230E-6;
tp = 5.4E-9;
o1 = 1E4;
o2 = 2E4;
o3 = 3E4;
o4 = 4E4;
o5 = 11E4;
k1 = 5.4E-3;
k2 = 10.8E-3;
k3 = 16.2E-3;
k4 = 22E-3;
k5 = 54E-3;
dy(1) = (P - y(1).*((abs(y(6)))^2 +1))./tc;
dy(2) = (P - y(2).*((abs(y(7)))^2 +1))./tc;
dy(3) = (P - y(3).*((abs(y(8)))^2 +1))./tc;
dy(4) = (P - y(4).*((abs(y(9)))^2 +1))./tc;
dy(5) = (P - y(5).*((abs(y(10)))^2 +1))./tc;
dy(6)= (y(1)-a).*((y(6))./tp) + (k1./tp).*(y(7)).*cos(y(11)) + (k5./tp).*(y(10))*cos(y(15));
dy(7)= (y(2)-a).*((y(7))./tp) + (k2./tp).*(y(8)).*cos(y(12)) + (k1./tp).*(y(6))*cos(y(11));
dy(8)= (y(3)-a).*((y(8))./tp) + (k3./tp).*(y(9)).*cos(y(13)) + (k2./tp).*(y(7))*cos(y(12));
dy(9)= (y(4)-a).*((y(9))./tp) + (k4./tp).*(y(10)).*cos(y(14)) + (k3./tp).*(y(8))*cos(y(13));
dy(10)= (y(5)-a).*((y(10))./tp) + (k5./tp).*(y(9)).*cos(y(14)) +(k4./tp).* (y(6))*cos(y(15));
dy(11) = o1 - (k1./tp).*(y(6)/y(7)).*(sin(y(11))) - (k1./tp).*(y(10)/y(6)).*(sin(y(11))) + (k2./tp).*(y(8)./y(7))*sin(y(12)) + (k5./tp).*(y(10)/y(6)).*sin(y(15));
dy(12) = o2 - (k2./tp).*(y(7)/y(8)).*(sin(y(12))) - (k2./tp).*(y(8)/y(7)).*(sin(y(12))) + (k3./tp).*(y(9)./y(8))*sin(y(13)) + (k1./tp).*(y(6)/y(7)).*sin(y(11));
dy(13) = o3 - (k3./tp).*(y(8)/y(9)).*(sin(y(13))) - (k3./tp).*(y(9)/y(8)).*(sin(y(13))) + (k4./tp).*(y(10)./y(9))*sin(y(14)) + (k2./tp).*(y(7)/y(8)).*sin(y(12));
dy(14) = o4 - (k4./tp).*(y(9)/y(10)).*(sin(y(14))) - (k4./tp).*(y(10)/y(9)).*(sin(y(14))) + (k5./tp).*(y(6)./y(10))*sin(y(15)) + (k3./tp).*(y(8)/y(9)).*sin(y(13));
dy(15) = o5 - (k5./tp).*(y(10)/y(6)).*(sin(y(15))) - (k5./tp).*(y(6)/y(10)).*(sin(y(15))) + (k1./tp).*(y(7)./y(6))*sin(y(11)) + (k4./tp).*(y(9)/y(10)).*sin(y(14));
end
0 Kommentare
Akzeptierte Antwort
Sam Chak
am 22 Aug. 2022
Hi @SAHIL SAHOO
Your k vector is assigned incorrectly. Try if this code works for you.
ti = 0;
tf = 1E-3;
tspan = [ti tf];
k = linspace(1, 2*pi, 5); % create 5 evenly-spaced points between 1 and 2pi
for i = 1:length(k)
y0 = [1; 1; 1; 1; 1; 1; 1; 1; 1; 1; k(i); 0; 0; 0; 0]; % just need y0 for each iteration for passing to ode45
tp = 5.4E-9;
[T,Y] = ode45(@(t,y) rate_eq(t, y), tspan, y0);
E(i) = mean(Y(:,11));
q(i) = (5.*(log(E(i))))./(2*pi);
plot(T, Y(:,11)), hold on % just to test if there are multiple plots
end
hold off
function dy = rate_eq(t, y)
dy = zeros(3*5, 1);
% parameters
P = 1.67;
a = 0.1;
tc = 230E-6;
tp = 5.4E-9;
o1 = 1E4;
o2 = 2E4;
o3 = 3E4;
o4 = 4E4;
o5 = 11E4;
k1 = 5.4E-3;
k2 = 10.8E-3;
k3 = 16.2E-3;
k4 = 22E-3;
k5 = 54E-3;
% dynamics
dy(1) = (P - y(1).*((abs(y(6)))^2 +1))./tc;
dy(2) = (P - y(2).*((abs(y(7)))^2 +1))./tc;
dy(3) = (P - y(3).*((abs(y(8)))^2 +1))./tc;
dy(4) = (P - y(4).*((abs(y(9)))^2 +1))./tc;
dy(5) = (P - y(5).*((abs(y(10)))^2 +1))./tc;
dy(6) = (y(1)-a).*((y(6))./tp) + (k1./tp).*(y(7)).*cos(y(11)) + (k5./tp).*(y(10))*cos(y(15));
dy(7) = (y(2)-a).*((y(7))./tp) + (k2./tp).*(y(8)).*cos(y(12)) + (k1./tp).*(y(6))*cos(y(11));
dy(8) = (y(3)-a).*((y(8))./tp) + (k3./tp).*(y(9)).*cos(y(13)) + (k2./tp).*(y(7))*cos(y(12));
dy(9) = (y(4)-a).*((y(9))./tp) + (k4./tp).*(y(10)).*cos(y(14)) + (k3./tp).*(y(8))*cos(y(13));
dy(10) = (y(5)-a).*((y(10))./tp) + (k5./tp).*(y(9)).*cos(y(14)) +(k4./tp).* (y(6))*cos(y(15));
dy(11) = o1 - (k1./tp).*(y(6)/y(7)).*(sin(y(11))) - (k1./tp).*(y(10)/y(6)).*(sin(y(11))) + (k2./tp).*(y(8)./y(7))*sin(y(12)) + (k5./tp).*(y(10)/y(6)).*sin(y(15));
dy(12) = o2 - (k2./tp).*(y(7)/y(8)).*(sin(y(12))) - (k2./tp).*(y(8)/y(7)).*(sin(y(12))) + (k3./tp).*(y(9)./y(8))*sin(y(13)) + (k1./tp).*(y(6)/y(7)).*sin(y(11));
dy(13) = o3 - (k3./tp).*(y(8)/y(9)).*(sin(y(13))) - (k3./tp).*(y(9)/y(8)).*(sin(y(13))) + (k4./tp).*(y(10)./y(9))*sin(y(14)) + (k2./tp).*(y(7)/y(8)).*sin(y(12));
dy(14) = o4 - (k4./tp).*(y(9)/y(10)).*(sin(y(14))) - (k4./tp).*(y(10)/y(9)).*(sin(y(14))) + (k5./tp).*(y(6)./y(10))*sin(y(15)) + (k3./tp).*(y(8)/y(9)).*sin(y(13));
dy(15) = o5 - (k5./tp).*(y(10)/y(6)).*(sin(y(15))) - (k5./tp).*(y(6)/y(10)).*(sin(y(15))) + (k1./tp).*(y(7)./y(6))*sin(y(11)) + (k4./tp).*(y(9)/y(10)).*sin(y(14));
end
0 Kommentare
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Historical Contests 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!