Plot the Bifurcation graph for a model equation.

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

Answer 1

[Vishnu]

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.

[ELISHA ANEBI]

Wow, thanks alot.

Answer 2

[khalid]

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')

[ELISHA ANEBI]

Thank you Khalid,

I am still trying to compute Bifurcation parameters from my ODE system of equation. If you can help my email is [e-mail]. I will appreciate any material that can help too