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

3 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

Categories

Find more on MATLAB in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!