Part A
 

(1) Write Matlab codes to generate the logistic curve by solving the logistic differential equation numerically (Please refer to the lecture notes on how to solve ordinary differential equations numerically). Compare the result with the analytical solution given in the lecture note. The logistic differential equation is given by

.

You may use a=b=1.0 and try two different initial condition P(0)=0.5 and P(0)=2.5. The comparison can be made by plotting the analytical result and the result obtained by solving the differential equation directly.

Ans.:

Create the M-files dgrowth1.m:

function dP = dgrowth1(t,P)
dP = (1.0-1.0*P)*P;

The following Matlab codes solve the equations with initial conditions P(0)=0.5 and P(0)=2.5 and compare with analytical results

% solving with a numerical method (ode45)
P0 = 0.5;
tspan = 0:0.1:5;
[t, P] = ode45('dgrowth1', tspan, P0);
plot(t, P, 'r-')
hold on

% compare with analytical result
Pa=P0*exp(t)./(1.0+P0*(exp(t)-1));
plot(t,Pa,'ro')
hold on

P0=2.5; %another initial condition
[t, P]= ode45('dgrowth1', tspan, P0);
plot(t, P, 'b-')
hold on

%compare with analytical result
Pa=P0*exp(t)./(1.0+P0*(exp(t)-1));
plot(t,Pa,'bo')



The numerical solutions of the ODE agree with the corresponding analytical results. The graph also shows that the solutions approach a fixed equilibrium value irrespective of the initial conditions.

(2) Consider the logistic map. As you increase the parameter r, you will see a sequence of bifurcation from a fixed point behavior to 2-cycle, then to 4-cycle, 8-cycle, etc. Find the values of the parameter r, where you see the bifurcation from a fixed point to a 2-cycle and from a 2-cycle to a 4-cycle. You may use the following codes to simulate the logistic map

r=2.5   %control parameter, r
x=0.7;  %the initial x is assigned to 0.7
xold=x; %initialize xold 

for k=1:1000 %iterate the map 1000 times
  xnew=r*xold*(1-xold); %new value
  x=[x xnew];           %add the new values to vector x
  xold=xnew;            %the new value becomes old value in the next step
end

Ans.:
We can go through the r value systematically (here we systematically change values of r from 2.0 to 3.5). By plotting or printing out the values of x in the last few iterations, we can determine where bifurcation occurs based on these values of x. The following program plots out the bifurcation diagram.

for r=2.5:0.05:3.5 %varies control parameter r systematically
  x=0.7; %the initial x is assigned to 0.7
  xold=x;

  for k=1:1000            %iterate the map 1000 times
    xnew=r*xold*(1-xold); %new value
    x=[x xnew];           %add the new values to vector x
    xold=xnew;            %the new value becomes old value in the next step
  end

  plot(r*ones(10), x(991:1000), 'bo') %plot values of x in the last 10 iterations for a given r.
  hold on
end

hold off



What are plotted in the above figure are values of x in the last 10 iterations (out of total 1000 iterations) for a given value of r. For a fixed point, there is only one x value; for the case of a 2-cycle, there are two distinct values for x; and for the case of a 4-cycle, x cycles through four distinct values. From the figure, it is easy to see that the first bifurcation (from fixed point to 2-cycle) occurs at r around 3.00, while the second bifurcation (from 2-cycle to 4-cycle) occurs at r around 3.45.

(3) Consider the case r=4 for the logistic map, generate a histogram for the probability distribution of x values.
 

Ans.:

r=4.0   %control parameter, r
x=0.7;  %the initial x is assigned to 0.7
xold=x; %initialize xold 

for k=1:1000 %iterate the map 1000 times
  xnew=r*xold*(1-xold); %new value
  x=[x xnew];           %add the new values to vector x
  xold=xnew;            %the new value becomes old value in the next step
end

hist(x)



(4) Generate the fern using the M-file (Iterated function systems) listed in the lecture note. What do you get if you change some parameters used in IFS slightly? (Give an example). Note that due to large number of iterations used, the it takes quite some time to generate a figure.
 

Ans.:  From the program listed in the lecture note,  we get
 

You can change the parameters in IFS and get a different fern. Here is the program that generates a slightly different fern.
For those who are curious: Note that there is a slight change from the program listed in the lecture note; the version of the program listed below runs faster because the memory for saving the matrix pts is allocated once from the beginning). A smaller MarkerSize is also used.

m1 = [0.85 0.14; -0.04 0.85];
c1 = [0 1.6]';
m2 = [-0.15 0.28; 0.26 0.24];
c2 = [0 0.44]';
m3 = [0.2 -0.26; 0.23, 0.22];
c3 = [0 1.6]'; m4 = [0 0; 0 0.16];
c4 = [0 0]'; x=[0 0]';
pts = zeros(2,10000);
for k = 1:10000
  t=rand(1,1);
  if t<=0.85 x=m1*x+c1;
  elseif t<=0.92 x=m2*x+c2;
  elseif t<=0.99 x=m3*x+c3;
  else x=m4*x+c4;
end
pts(:,k) = x;
end
plot(pts(1,:),pts(2,:),'.', 'MarkerSize',2)
axis equal

Part B (optional) One-Dimensional Random walk (this might be difficult) Generate 100 one-dimensional random walks: x(n+1)=x(n)+d(n), where d(n) is a random number which takes values 1 and -1 with equal probability. Now evaluate <x(n)> and <(x(n))²>, where the averages (indicated by <...>) are over all 100 walks. How do <x(n)> and <(x(n))²> depend on n.

Ans.:

N=1000;
p=2*rand(N,100)-1;
walk=cumsum([zeros(1,100);p]);
a=mean(walk');
b=mean((walk.*walk)');
plot(a)
hold on
plot(b)

One can see that <(x(n))²> is roughly proportional to n and <x(n)> is close to zero.