Demo of MATLAB using the example of Bifurcation

From Math Images
Jump to: navigation, search

Problem Statement

Known:  x_{n+1} = r(1-x_n)x_n where  x_n is a number between  0 and  1 .

The Claim: we can plot the graph as follows regardless of  x_0 value

LogisticMap BifurcationDiagram.png

where the x axis is the  r value and the  y axis is the value of  x_n when  n is very big.

Solution Numerical Simulation: We want to write a program that takes an initial value,  x_0 , and compute  \lim_{n \to \infty}x_n for every  r value between  0 and some certain value.

Trial Test: Let's try something small first with MATLAB and see if the values computed agree with the picture at all.

MATLAB Code Explanation
 %   let's clear the command window first

 %   declare constant first
 %   variable names have to be meaningful
 x_0 = 0.01;
 r = 2.8;

 %   create an empty vector of size 1x1000
 x_vec = zeros(1,1000);

 %   let the first element of the vector be
 %   x_0
 x_vec(1) = x_0;

 %   now all we have to do is to iterate 1000 times
 for i = 2:1000
     x_vec(i) = r*x_vec(i-1)*(1-x_vec(i-1));

 %   let's display the last element (1000th) of the vector
 %   and that is the x_n value when n -> inf
 disp(sprintf('Taking lim n -> inf, r = %f and x_0 = %f, we have x_n = %f',r, x_0, x_vec(1000)));

and the output of the program is

Taking lim n -> inf, r = 2.800000 and x_0 = 0.010000, we have x_n = 0.642857

Anything that follows a % in MATLAB script is ignored by the program. Thus all comments about the program follows the % sign. It reminds what your intentions were when you wrote the program so that when you come back to it some time later, you can still understand it.


Each line of code ends with ;. If you don't end it with ;, MATLAB will output the results to the command window. Sometimes, it is what you intended. However, you want to suppress the output until all calculations are finished.


Clears command window. Basically it clears the output window.

x_0 = 0.01

All variables have to have meaningful names, to you. You don't have to specify what type of the data the variable is like you would do in C/C++. MATLAB figures it out itself. So a variable can be an integer, float (i.e. 1.2355), string (i.e. MATLAB), complex number (i.e. 5i) and etc.

x_vec = zeros(1,1000)

This creates a 1x1000 vector of integer zeros. The purpose of this is to fast computing. MATLAB does not have to dynamically allocate new memory space every time we have a new x_i value. It just have to erase the zero and replace it with the computed value.

x_vec(1) = x_0;

This is how we index into a vector. The first element of x_vec is x_vec(1) and we want it to equal to x_0.

for i = 2:1000
     x_vec(i) = r*x_vec(i-1)*(1-x_vec(i-1));

This is the for loop. This means that the value of i starts from 2 and increments to 1000 by step of 1. Everytime the loop iterates, the code in the middle executes. What we are doing here is simply updating the x_n value and stop at x_{1000}.

Don't worry about the last line of the code. It is just to output the value to the command window.

Now this program works for r=2.8. What about r=3.2 where there is two x_\infty values? Well, we can simply change the r value and run the simulation again. We then get

Taking lim n -> inf, r = 3.200000 and x_0 = 0.010000, we have x_n = 0.799455

This certainly agrees with the upper part of the graph. What about the lower one? Now we need to give it a new  x_0 value because with the same  x_0 value, now matter how many times you run the program, you will get the same result. So with a new  x_0 value, we have

Taking lim n -> inf, r = 3.200000 and x_0 = 0.850000, we have x_n = 0.513045

This value agrees with the picture. So that means our code is working fine. Now all we have to do is to iterate this program for every value of r and lots of different  x_0 values. Eventually we can then plot the picture with all those values.

Full Code

%   creates a 1x1000 vector of values evenly spaces between 0 and 3.6 
r_vec = linspace(0,3.6,1000);
%   creats an empty 1x1000 vector to hold the values of every iteration
x_vec = zeros(1,1000);
%   creates an empty 1x32 vector to hold the values of every initial
%   condistions. by right we need a lot more since when r values gets
%   really big, x values changes between exponentially increasing amount of
%   numbers but we are stopping at r = 3.6 so we don't have to worry about
%   that
x_temp = zeros(1,32);

%   for every r value
for j = 1:1000        
    %   try different initial conditions
    for k = 1:32
        x_vec(1) = rand;
        %   and interate 1000 times to get long term values
        for i=1:999
            x_vec(i+1) = r_vec(j)*x_vec(i)*(1-x_vec(i));
        x_temp(k) = x_vec(1000);
        hold on;

The code presented above is neither run-time efficient nor memory space efficient. But it gets the job done. Here is the picture.

OneK rvalue.png

I will not explain this code in detail like I did previously. If you guys are interested, you or we can work together to write a more efficient program to plot over larger range of r values.