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

2 views (last 30 days)
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 Comment
Jan
Jan on 8 Jul 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.

Sign in to comment.

Accepted Answer

KSSV
KSSV on 8 Jul 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

More Answers (1)

Hrusheekesh
Hrusheekesh on 11 Jul 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

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!