Demo of MATLAB using the example of Bifurcation
Known: where is a number between and .
The Claim: we can plot the graph as follows regardless of value
where the axis is the value and the axis is the value of when is very big.
Solution Numerical Simulation: We want to write a program that takes an initial value, , and compute for every value between 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.
% let's clear the command window first clc; % 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)); end % 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. ), string (i.e. ), complex number (i.e. ) 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 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)); end
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 value and stop at .
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 . What about where there is two values? Well, we can simply change the 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 value because with the same value, now matter how many times you run the program, you will get the same result. So with a new 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 and lots of different values. Eventually we can then plot the picture with all those values.
clc; clf; % 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)); end x_temp(k) = x_vec(1000); hold on; plot(r_vec(j),x_temp(k)); end end
The code presented above is neither run-time efficient nor memory space efficient. But it gets the job done. Here is the picture.
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 values.