if elseif else statement not working
1 view (last 30 days)
Show older comments
I have this code
r0=0.05;
k1=0.5;
k2=0.5;
K=1;
mu=0.5;
rho=0.5;
alpha=0.1;
b=0.05;
epsilon=0.25;
sigma=0.05;
eta=0.05;
pi=0.05;
omega0=0.001;
syms M m Msol positive
P(M) = epsilon*M/(1+M);
N(M) = mu+(rho*M/(1+M));
lambda(M)=1/(1+N);
r(M) = r0*(1+k1*N*(1-k2*N));
F = (m/b)*r(M)*(eta+P(M))*(1+N(M))*(1-(((alpha*sigma*m*(eta+P(M))*(1+N(M))+b*(m+alpha)*(pi*M-omega0)))/(b*alpha*sigma*K)))-(m/sigma)*(pi*M-omega0)==0;
[num,den]=numden(lhs(F));
m_num = linspace(0.01,0.9,50);
for u = 1:numel(m_num)
Msol(u) = vpa(solve(subs(num,m,m_num(u))));
end
R0 = (b.*K.*lambda(omega0./pi))./(m_num.*(eta+P(omega0./pi)));
if R0 < 1
l = 0;
y = 0;
v = 0;
else
l = ((b.*K.*lambda(omega0./pi).*(pi.*Msol-omega0))./(alpha.*sigma.*R0.*(eta+P(omega0./pi))));
y = ((pi.*Msol-((pi+omega0)./epsilon).*(((b.*K.*lambda(omega0./pi))./(m_num.*R0))-eta))./(sigma));
v = ((m_num.*R0.*(eta+P(omega0./pi)).*(pi.*Msol-omega0))./(sigma.*K.*lambda(omega0./pi).*(eta+P(Msol))));
end
subplot(3,1,1)
p1 = plot (R0,l);
ylabel('l^*')
xlim([0 5]);
ylim([-0.15 0.15]);
grid on;
subplot(3,1,2)
p2 = plot (R0,y);
ylabel('y^*')
xlim([0 5]);
ylim([-0.15 0.15]);
grid on;
subplot(3,1,3)
p3 = plot (R0,v);
xlabel('R_0'), ylabel('v^*')
xlim([0 5]);
ylim([-0.15 0.15]);
grid on;
set(p1,{'LineWidth'},{1.5});
set(p2,{'LineWidth'},{1.5});
set(p3,{'LineWidth'},{1.5});
I am trying to use if,elseif,else syntax to have l,y and v zeroes for R0<1 and positive elsewhere. However, it looks like MATLAB does not recognize this part. Any help would be appreciated!
0 Comments
Accepted Answer
Star Strider
on 25 Apr 2023
Edited: Star Strider
on 25 Apr 2023
Use a logical vector to select the elements you want to plot, then use that as a subscript to both arguments in each plot call —
r0=0.05;
k1=0.5;
k2=0.5;
K=1;
mu=0.5;
rho=0.5;
alpha=0.1;
b=0.05;
epsilon=0.25;
sigma=0.05;
eta=0.05;
pi=0.05;
omega0=0.001;
syms M m Msol positive
P(M) = epsilon*M/(1+M);
N(M) = mu+(rho*M/(1+M));
lambda(M)=1/(1+N);
r(M) = r0*(1+k1*N*(1-k2*N));
F = (m/b)*r(M)*(eta+P(M))*(1+N(M))*(1-(((alpha*sigma*m*(eta+P(M))*(1+N(M))+b*(m+alpha)*(pi*M-omega0)))/(b*alpha*sigma*K)))-(m/sigma)*(pi*M-omega0)==0;
[num,den]=numden(lhs(F));
m_num = linspace(0.01,0.9,50);
for u = 1:numel(m_num)
Msol(u) = vpa(solve(subs(num,m,m_num(u))));
end
R0 = (b.*K.*lambda(omega0./pi))./(m_num.*(eta+P(omega0./pi)));
l = ((b.*K.*lambda(omega0./pi).*(pi.*Msol-omega0))./(alpha.*sigma.*R0.*(eta+P(omega0./pi))));
y = ((pi.*Msol-((pi+omega0)./epsilon).*(((b.*K.*lambda(omega0./pi))./(m_num.*R0))-eta))./(sigma));
v = ((m_num.*R0.*(eta+P(omega0./pi)).*(pi.*Msol-omega0))./(sigma.*K.*lambda(omega0./pi).*(eta+P(Msol))));
Lv = R0 >= 1; % Logical Vector
subplot(3,1,1)
p1 = plot (R0(Lv),l(Lv));
hold on
plot([0 1.02], [0 0], '-', 'Color',p1.Color, 'LineWidth',1)
hold off
ylabel('l^*')
xlim([0 5]);
ylim([-0.15 0.15]);
grid on;
subplot(3,1,2)
p2 = plot (R0(Lv),y(Lv));
hold on
plot([0 1.02], [0 0], '-', 'Color',p2.Color, 'LineWidth',1)
hold off
ylabel('y^*')
xlim([0 5]);
ylim([-0.15 0.15]);
grid on;
subplot(3,1,3)
p3 = plot (R0(Lv),v(Lv));
hold on
plot([0 1.02], [0 0], '-', 'Color',p3.Color, 'LineWidth',1)
hold off
xlabel('R_0'), ylabel('v^*')
xlim([0 5]);
ylim([-0.15 0.15]);
grid on;
set(p1,{'LineWidth'},{1.5});
set(p2,{'LineWidth'},{1.5});
set(p3,{'LineWidth'},{1.5});
EDIT — (25 Apr 2023 at 16:54)
Added constant lines from [0 1.02] in every plot.
.
4 Comments
More Answers (0)
See Also
Categories
Find more on Axis Labels 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!