I want to need (maxY(:,1)) for different value of k, how the for loop will help?

2 vues (au cours des 30 derniers jours)
I need the "max(Y(:,1))" for the different value of k, in this code I put a single constant value of k and I need the max(Y(:,1) for the inteval of k =[0:0.1E-3 : 3E-3], I'm not able to get an idea of putting a for loop to get this.
ti = 0;
tf = 100E-4;
tspan=[ti tf];
y0=[1;1;1;1;1].*10E-2;
[T,Y]= ode45(@(t,y) rate_eq(t,y),tspan,y0);
subplot 211
plot(T,Y(:,2));
ylim([0 5])
xlim([0 0.3E-3])
subplot 212
plot(T,((Y(:,5))./6.28));
xlim([0 0.3E-3])
disp(max(Y(:,1)))
function dy = rate_eq(t,y)
dy = zeros(5,1);
o = 2E5;
tc = 30E-9;
tf = 230E-6;
a1 = 0.1;
a2 = 0.1;
P1 = 0.2;
P2 = 0.2;
k= 2E-3
V = 1;
y11 = y(1);
y22 = y(2);
y33 = y(3);
y44 = y(4);
y55 = y(5);
dy(1) = (1/tc).*((y33)- a1).*(y(1)./V) +(k./tc).*(y(2)./V).*cos((y55));
dy(2) = (1/tc).*((y44)- a2).*(y(2)./V) +(k./tc).*((y11)./V).*cos((y55));
dy(3) = (1/tf).*(P1 - (y33).*(abs(y11).^2 +1));
dy(4) = (1/tf).*(P2 - (y44).*(abs(y22).^2+1));
dy(5) = o - (k./tc).*(((y11)./(y22)) + ((y22)./(y11)).*sin((y55)));
end
  1 commentaire
Jan
Jan le 8 Juil 2022
Some hints:
subplot 211
This notation is deprecated for 20 years now. Prefer: subplot(2,1,1) or use the modern tiledlayout.
abs(x).^2 is the same as x^2 for real scalar x. There is no need to include a single variable in parentheses. The alis-variables y11, y22... are not useful here. Compare:
y11 = y(1);
y22 = y(2);
y33 = y(3);
y44 = y(4);
y55 = y(5);
dy(1) = (1/tc).*((y33)- a1).*(y(1)./V) +(k./tc).*(y(2)./V).*cos((y55));
dy(2) = (1/tc).*((y44)- a2).*(y(2)./V) +(k./tc).*((y11)./V).*cos((y55));
dy(3) = (1/tf).*(P1 - (y33).*(abs(y11).^2 +1));
dy(4) = (1/tf).*(P2 - (y44).*(abs(y22).^2+1));
dy(5) = o - (k./tc).*(((y11)./(y22)) + ((y22)./(y11)).*sin((y55)));
with the readability of:
dy(1) = ((y(3) - a1) * y(1) / V + k * y(2) / V * cos(y(5))) / tc;
dy(2) = ((y(4) - a2) * y(2) / V + k * y(1) / V * cos(y(5))) / tc;
dy(3) = (P1 - y(3) * (y(1)^2 + 1)) / tf;
dy(4) = (P2 - y(4) * (y(2)^2 + 1)) / tf;
dy(5) = o - k / tc * (y(1) / y(2) + y(2) / y(1) * sin(y(5)));
The clearer the code, the easier you see typos.

Connectez-vous pour commenter.

Réponse acceptée

KSSV
KSSV le 8 Juil 2022
ti = 0;
tf = 100E-4;
tspan=[ti tf];
y0=[1;1;1;1;1].*10E-2;
k =[0:0.1E-3 : 3E-3] ;
iwant = zeros(length(k),1) ;
for i = 1:length(k)
[T,Y]= ode45(@(t,y) rate_eq(y,k(i)),tspan,y0);
iwant(i) = max(Y(:,1)) ;
end
iwant
iwant = 31×1
2.3962 2.3941 2.3856 2.3666 2.3504 2.3302 2.2943 2.2955 2.3370 2.3632
function dy = rate_eq(y,k)
dy = zeros(5,1);
o = 2E5;
tc = 30E-9;
tf = 230E-6;
a1 = 0.1;
a2 = 0.1;
P1 = 0.2;
P2 = 0.2;
% k= 2E-3
V = 1;
y11 = y(1);
y22 = y(2);
y33 = y(3);
y44 = y(4);
y55 = y(5);
dy(1) = (1/tc).*((y33)- a1).*(y(1)./V) +(k./tc).*(y(2)./V).*cos((y55));
dy(2) = (1/tc).*((y44)- a2).*(y(2)./V) +(k./tc).*((y11)./V).*cos((y55));
dy(3) = (1/tf).*(P1 - (y33).*(abs(y11).^2 +1));
dy(4) = (1/tf).*(P2 - (y44).*(abs(y22).^2+1));
dy(5) = o - (k./tc).*(((y11)./(y22)) + ((y22)./(y11)).*sin((y55)));
end

Plus de réponses (1)

Hrusheekesh
Hrusheekesh le 11 Juil 2022
Hi sahil,
what you can do is send the values of k in a array.
k =[0:0.1E-3:3E-3];
r=[];
for i = 1:size(k)
[T,Y]= ode45(@(t,y) rate_eq(y,k(i)),tspan,y0);
r = [r max(Y(:,1))];
end

Catégories

En savoir plus sur Logical 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!

Translated by