Bifurcation Diagram

I want to draw the bifurcation diagram fro the model.

dx/dt=rx(1-x/K)-mxy/(ax+by+c) 
dy/dt=emxy/(ax+by+c)-dy-hy^2.

parameters are all +ve.

I have tryed to plot it but fails.

clear
r=0.806; a=15; b=16;c=17;e=0.333;d=0.3;h=0.01;K=200;
x (1)=0.7;
y (1)=0.11;
t (1)=0;
for m=6:1:22
for i=1:10000
    t (i+1)=t (i)+.01;
    x (i+1)= x (i)+0.01*[r*x (i)*(1-x (i)/K)-m*x (i)*y (i)/(a*x \
(i)+b*y (i)+c)];
   y (i+1)=y (i)+.01*[e*m*x (i)*y (i)/(a*x (i)+b*y (i)+c)-d*y (i)-h*y \
(i)^2];
  end
plot (m,x, 'b')
hold on; 
end
xlabel ('m')
ylabel ('x')
figure (2)
for m=6:1:22
for i=1:10000
    t (i+1)=t (i)+.01;
    x (i+1)= x (i)+0.01*[r*x (i)*(1-x (i)/K)-m*x (i)*y (i)/(a*x \
(i)+b*y (i)+c)];
   y (i+1)=y (i)+.01*[e*m*x (i)*y (i)/(a*x (i)+b*y (i)+c)-d*y (i)-h*y \
(i)^2];
  end
plot (m,y, 'b')
hold on; 
end
xlabel ('m')
ylabel ('y')

Please modify or help me to modify the matlab code to draw the following bifurcation diagram (parameter VS population):

1.Transcritical bifurcation (x vs m & y vs. m) around at m= 13.666

2. Saddle-node bifurcation (x vs m & y vs. m) around at m = 20.8.

3. Hopf-bifurcation (x vs m & y vs. m) at m=14.73, (d,h) = (0.02,0.001) and others are same.

[Mirza]

Plz try to read this. It has all relevant information. http://amath.colorado.edu/courses/2460/2009fall/Homeworks/bifurcation.pdf

[DEEPIHA PADMANATHAN]

The requested page "/courses/2460/2009fall/Homeworks/bifurcation.pdf" could not be found. Page not found

[Walter Roberson]

https://amath.colorado.edu/sites/default/files/2015/03/274176680/APPM2460Worksheet11.pdf

[MD EASHIN PATWARY]

Both file is not working,

Answer 1

[Rick Rosson]

Hi Pallav,

There are several issues with the code that you posted. Some of these problems are syntax errors, some are non-standard coding style, and some are non-optimal formatting.

Syntax Errors

The biggest syntax error, and the main reason this code fails to run, is the improper use of square brackets [ and ] in places where parentheses (| and |) are required.

So, as a first step, please replace every instance of [ with (|, and likewise, every instance of |] with ).

MATLAB does use square brackets for certain purposes, but never as a way of grouping expressions within assignment statements to specify the order of evaluation. MATLAB always uses parentheses for that purpose, never square brackets.

Formatting

Please do not put a blank space between the name of an array variable and the open parenthesis that precedes the index of the array. Also, you should include at least one blank space before and after the = sign in any assignment statement. Likewise, you should also include at least one blank space on either side of any + or - sign in an expression consisting of more than one term.

So, for example, instead of writing

   t (i+1)=t (i)+.01;

please write

   t(i+1) = t(i) + 0.01;

Although the first line is perfectly valid MATLAB syntax, it it horribly difficult to read and understand, whereas the second line is far easier to read and understand. The first line is also generally considered contrary to most MATLAB coding standards and conventions, whereas the second line is generally considered consistent with standard conventions.

Parameterization

It is often helpful to introduce one or more parameters into the code as a way of improving the readability and maintainability of the code. So, for example, instead of using the literal number 0.01 at several places throughout the code, it may make sense to define a parameter dt (meaning "delta-t" or "time increment") at the beginning of the script:

   dt = 0.01;

and then replacing every instance of 0.01 with dt, such as:

   t(i+1) = t(i) + dt;

and

   x(i+1) = x(i) + dt*( ... );

and so on.

Pre-Allocation

It is very important to pre-allocate time t and the two state variables x and y to avoid having them grow during each iteration of the inner for loop. To do so, please insert the following four lines of code after the outer for loop and prior to the inner loop:

    N = 10000;
    x = zeros(N,1);
    y = zeros(N,1);
    t = zeros(N,1);

Although this code is not strictly necessary, it will reduce the run-time of your script quite substantially, and it will avoid unnecessarily fragmenting the memory during execution of the script.

Initial Conditions

Now that you have pre-allocated the time and state variables, you need to set the initial conditions for each of these three variables. So please cut the following three lines of code and paste them after the pre-allocation code and prior to the inner for loop:

    x(1) = 0.7;
    y(1) = 0.11;
    t(1) = 0;

By moving these three lines inside the outer for loop but prior to the inner for loop, you will ensure that each instance of the parameter m will start with the correct initial conditions independent of the prior instance.

Modified Code

Here is what the code looks like after all of the changes mentioned above (plus a few others):

    %%Initialize the environment

    close all;
    clear all;
    clc;

    %%Model parameters

    r = 0.806; 
    a = 15; 
    b = 16;
    c = 17;
    e = 0.333;
    d = 0.3;
    h = 0.01;
    K = 200;

    %%Time parameters

    dt = 0.01;
    N = 10000;

    %%Set-up figure and axes

    figure;

    ax(1) = subplot(2,1,1);
    hold on
    xlabel ('m');
    ylabel ('x');

    ax(2) = subplot(2,1,2);
    hold on
    xlabel ('m');
    ylabel ('y');

    %%Main loop

    for m=6:1:22

        x = zeros(N,1);
        y = zeros(N,1);
        t = zeros(N,1);

        x(1) = 0.7;
        y(1) = 0.11;
        t(1) = 0;

        for i=1:N
            t(i+1) = t(i) + dt;
            x(i+1) = x(i) + dt*(  r*x(i)*(1-x(i)/K)-m*x(i)*y(i)/(a*x (i)+b*y(i)+c)  );
            y(i+1) = y(i) + dt*(  e*m*x(i)*y(i)/(a*x(i)+b*y(i)+c)-d*y(i)-h*y(i)^2  );
        end

        plot(ax(1),m,x,'color','blue','marker','.');
        plot(ax(2),m,y,'color','blue','marker','.');

    end

HTH.

Rick

[Atom]

Thank you for the suggestions. I have learned and changed the bracket issue according to your suggestion. It has some more error as you have mentioned, and for this reason the code is still not working. Hope you we suggest further to improve the code at your earliest convenience.

[Atom]

Thank you again. I have corrected the =,-,+ sign issues according to your suggestions.

[Atom]

Thank you. I have written dt = 0.01; at the beginning of the script and all 0.01 has been replaced by dt .

[Atom]

Thank you very much for your detailed suggestions. It helps me to learn to make further codes. But, I am still unable to get the bifurcation diagrams. The fig that I am getting is not the bif diagram.

[Rizwana Junaid]

ur system equations are differential equations not the discrete equations. so u should try to solve ur equations with solver ode45. so introduce solver command in the m loop.