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??
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
SAHIL SAHOO
el 22 de Ag. de 2022
Respondida: Sam Chak
el 22 de Ag. de 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 comentarios
Respuesta aceptada
Sam Chak
el 22 de Ag. de 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 comentarios
Más respuestas (0)
Ver también
Categorías
Más información sobre Historical Contests en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!