Plot the Bifurcation graph for a model equation.
15 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
ELISHA ANEBI
am 11 Jul. 2023
Kommentiert: ELISHA ANEBI
am 16 Mai 2024
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 Kommentare
Akzeptierte Antwort
Vishnu
am 12 Jul. 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.
Weitere Antworten (1)
khalid
am 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')
Siehe auch
Kategorien
Mehr zu Nonlinear Dynamics finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!