Plot the Bifurcation graph for a model equation.
14 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Kindly help me bifurcation diagram for the equation and parameter value below. I have tried getting the graph but it is giving me graph out of range.
% The Matlab Codes for the forward bifurcation diagram
Rev_value=0.018:0.01:4;
Root_array=zeros (length (Rev_value), 2);
qI=0.001923; qA=0.00000004013; etaA=0.1213; etaQ=0.003808; w=0.5925;
Lambda=0.1598643e-7; theta=0.022;
delta=0.125; mu=0.01119; pi=0.464360344; deltaQ=6.847e-4; beta=1.086e-1;
qE=1.8113e-4; rhoQ=0.0815; a=0.16255; k=0.15; v1=0.71; v2=0.29;alpha=0.57e-1; deltaI=0.00000000223; rhoA=0.1; rhoI=0.0666666;
c1=mu+v1; c2=1-w; c3=mu+alpha+v2;c4=1-k; c5=qE+delta+mu; c6=rhoA+mu; c7=delta*(1-a); c8=rhoI+qI+deltaI+mu; c9=rhoA+deltaQ+mu;
b1=1-theta;
B=delta*a*c8*c9+c7*k*qI*rhoQ+c8*k*qE*rhoQ;
G=qI*c7*c6+qE*c8*c6;
H1=C2*c5*c6*c8*c9;
H2=c5*c6*c8*c9*(c3*+c1*c2)-(b1*c2+c2*theta)*(c6*c9*c7*pi*beta+pi*G*beta*etaQ+pi*B*beta*etaA);
H3=c5*c6*c8*c9*(c3*c1-v1*v2)*(1-b1*Rev_value)-(c3*theta+c2*v1*theta)*(c6*c9*c7*pi*beta+pi*G*beta*etaQ+pi*B*beta*etaA);
hold on
for i=1:1:length(Rev_value);
Rev=Rev_value(i);
% bifurcation parameter
%Coefficients of quadratic equation H1, H2, H3
p=[H1,H2,H3];
r =roots(p);
len=length(r);
for t=1:1:len
0 commentaires
Réponse acceptée
Vishnu
le 12 Juil 2023
Hi ELISHA ANEBI,
I noticed that in this equation 'H1=C2*c5*c6*c8*c9;' the C2 sholud be changed to c2. Additionally, there is a syntax error in the line for t=1:1:len;. The semicolon at the end should be removed.
To create the bifurcation diagram, you can modify the code as follows:
Rev_value=0.018:0.01:4;
Root_array=zeros (length (Rev_value), 2);
qI=0.001923; qA=0.00000004013; etaA=0.1213; etaQ=0.003808; w=0.5925;
Lambda=0.1598643e-7; theta=0.022;
delta=0.125; mu=0.01119; pi=0.464360344; deltaQ=6.847e-4; beta=1.086e-1;
qE=1.8113e-4; rhoQ=0.0815; a=0.16255; k=0.15; v1=0.71; v2=0.29;alpha=0.57e-1; deltaI=0.00000000223; rhoA=0.1; rhoI=0.0666666;
c1=mu+v1; c2=1-w; c3=mu+alpha+v2;c4=1-k; c5=qE+delta+mu; c6=rhoA+mu; c7=delta*(1-a); c8=rhoI+qI+deltaI+mu; c9=rhoA+deltaQ+mu;
b1=1-theta;
B=delta*a*c8*c9+c7*k*qI*rhoQ+c8*k*qE*rhoQ;
G=qI*c7*c6+qE*c8*c6;
H1=c2*c5*c6*c8*c9;
H2=c5*c6*c8*c9*(c3*+c1*c2)-(b1*c2+c2*theta)*(c6*c9*c7*pi*beta+pi*G*beta*etaQ+pi*B*beta*etaA);
H3=c5*c6*c8*c9*(c3*c1-v1*v2)*(1-b1*Rev_value)-(c3*theta+c2*v1*theta)*(c6*c9*c7*pi*beta+pi*G*beta*etaQ+pi*B*beta*etaA);
hold on
for i = 1:length(Rev_value)
Rev = Rev_value(i);
% Calculate G
G = qI * c7 * c6 + qE * c8 * c6;
% Coefficients of quadratic equation H1, H2, H3
H1 = c2 * c5 * c6 * c8 * c9;
H2 = c5 * c6 * c8 * c9 * (c3 + c1 * c2) - (b1 * c2 + c2 * theta) * (c6 * c9 * c7 * pi * beta + pi * G * beta * etaQ + pi * B * beta * etaA);
H3 = c5 * c6 * c8 * c9 * (c3 * c1 - v1 * v2) * (1 - b1 * Rev) - (c3 * theta + c2 * v1 * theta) * (c6 * c9 * c7 * pi * beta + pi * G * beta * etaQ + pi * B * beta * etaA);
% Coefficients of quadratic equation p = [H1, H2, H3]
p = [H1, H2, H3];
r = roots(p);
len = length(r);
for t = 1:len
% Store the real part of the roots in the Root_array
if isreal(r(t))
Root_array(i, t) = real(r(t));
end
end
end
% Plot the bifurcation diagram
plot(Rev_value, Root_array, '.')
xlabel('Rev')
ylabel('Roots')
title('Bifurcation Diagram')
The code will generate a bifurcation diagram by plotting the real part of the roots against the parameter Rev.
Plus de réponses (1)
khalid
le 15 Mai 2024
Rev_value=0.018:0.01:4;
Root_array=zeros (length (Rev_value), 2);
qI=0.001923; qA=0.00000004013; etaA=0.1213; etaQ=0.003808; w=0.5925;
Lambda=0.1598643e-7; theta=0.022;
delta=0.125; mu=0.01119; pi=0.464360344; deltaQ=6.847e-4; beta=1.086e-1;
qE=1.8113e-4; rhoQ=0.0815; a=0.16255; k=0.15; v1=0.71; v2=0.29;alpha=0.57e-1; deltaI=0.00000000223; rhoA=0.1; rhoI=0.0666666;
c1=mu+v1; c2=1-w; c3=mu+alpha+v2;c4=1-k; c5=qE+delta+mu; c6=rhoA+mu; c7=delta*(1-a); c8=rhoI+qI+deltaI+mu; c9=rhoA+deltaQ+mu;
b1=1-theta;
B=delta*a*c8*c9+c7*k*qI*rhoQ+c8*k*qE*rhoQ;
G=qI*c7*c6+qE*c8*c6;
H1=c2*c5*c6*c8*c9;
H2=c5*c6*c8*c9*(c3*+c1*c2)-(b1*c2+c2*theta)*(c6*c9*c7*pi*beta+pi*G*beta*etaQ+pi*B*beta*etaA);
H3=c5*c6*c8*c9*(c3*c1-v1*v2)*(1-b1*Rev_value)-(c3*theta+c2*v1*theta)*(c6*c9*c7*pi*beta+pi*G*beta*etaQ+pi*B*beta*etaA);
hold on
for i = 1:length(Rev_value)
Rev = Rev_value(i);
% Calculate G
G = qI * c7 * c6 + qE * c8 * c6;
% Coefficients of quadratic equation H1, H2, H3
H1 = c2 * c5 * c6 * c8 * c9;
H2 = c5 * c6 * c8 * c9 * (c3 + c1 * c2) - (b1 * c2 + c2 * theta) * (c6 * c9 * c7 * pi * beta + pi * G * beta * etaQ + pi * B * beta * etaA);
H3 = c5 * c6 * c8 * c9 * (c3 * c1 - v1 * v2) * (1 - b1 * Rev) - (c3 * theta + c2 * v1 * theta) * (c6 * c9 * c7 * pi * beta + pi * G * beta * etaQ + pi * B * beta * etaA);
% Coefficients of quadratic equation p = [H1, H2, H3]
p = [H1, H2, H3];
r = roots(p);
len = length(r);
for t = 1:len
% Store the real part of the roots in the Root_array
if isreal(r(t))
Root_array(i, t) = real(r(t));
end
end
end
% Plot the bifurcation diagram
plot(Rev_value, Root_array, '.')
xlabel('Rev')
ylabel('Roots')
title('Bifurcation Diagram')
Voir également
Catégories
En savoir plus sur Nonlinear Dynamics 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!