Biffurcation of the Lorenz Equation

Bifurcation is a way to draw the attractor of a chaotic system. You have some iterated system with a parameter, you vary the parameter and calculate the values the system takes. The goal is to plot the “stable values”, which is most easily achievable by simply iterating the sequence for a while and plotting a number of values that the system takes after it has been allowed to stabilize.

I decided that while this method may be fine for some systems, it isn’t capable of getting a decent picture in a less predictable system with a very strange attractor. For this we want a more probabilistic approach. Some sort of probability cloud which contains the system.

How much time does the system spend on any given area? I set out to create an image that would answer this question to one extent or another. The motivation for me to develop such imaging was the lorenz equation. What is commonly referred to as the Lorenz attractor, is the lorenz equation with , and . For these values the attractor is strange. This is one of the first examples of numerical chaos.

The lorenz equation is a coupled set of differential equations:

We can varry any of these constants, , , and plot the values x, y, or z takes against them. This is a basic explanation of how standard bifurcation is done. This is also often referred to as a Feigenbaund diagram.

But to get a more complete picture, we want to take into account how often the system visits any given value.

So we divide the domain of the variable in question, x, y, or z up into discreet intervals. Then we count the number of values that fall in each bin. Just like a histogram. We tally these values up in a matrix for plotting. Every value is color coded. I use the MATLAB imagesc() command, which takes a matrix and plots it as an image, higher values on the matrix corresponding to hotter colors. I had to filter the image with a logarithm to make it more clear, without suck filtering it was very faint.

We plot the x values vs where x is the horizontal axis.

Here it is for

i apologize, the axes are not labeled correctly. Here is a zoom in on the interface between the normal attractor and advent of the strange attractor.

More interesting images to come. Here is the code i used to generate these images:

m = 2000; % rho resolution
M = 2000; % X resolution
tf = 20; %final time
R = linspace(1,60,m);

for k=1:m
clear x y z;
x(1) = 3;
y(1) = 15;
z(1) = 5;
dt = 10^-2;
pran = 10;
ray = R(k);
phy = 8/3;
N = round(tf/dt);
dt = tf/N;
for j=1:N
    x1 = x(j) + dt*(pran*(y(j)-x(j)));
    x1 = 0.75*x(j) + 0.25*(x1 + dt*(pran*(y(j)-x1)));
        x(j+1) = (x(j) + 2*(x1 + dt*(pran*(y(j)-x1))))/3;

y1 = y(j) + dt*(x(j)*(ray-z(j))-y(j));
    y1 = 0.75*y(j) + 0.25*(y1 + dt*(x(j)*(ray-z(j))-y1));
        y(j+1) = (y(j) + 2*(y1 + dt*(x(j)*(ray-z(j))-y1) ))/3;

z1 = z(j) + dt*(x(j)*y(j) - phy*z(j));
    z1 = 0.75*z(j) + 0.25*(z1 + dt*(x(j)*y(j) - phy*z1));
        z(j+1) = (z(j) + 2*(z1 + dt*(x(j)*y(j) - phy*z1)))/3;

end
k/(2*m)
Z(k,:) = z;
Y(k,:) = y;
X(k,:) = x;
end

draw = zeros(M);

xmin = min(min(X));
xmax = max(max(X));

xx = linspace(xmax,xmin,M);
st = round(N/10);
for j=1:m
        for k=1:M-1
                draw(k,j) = sum((xx(k) >= X(j,st:end)).*(xx(k+1) <= X(j,st:end))
);
        end
(j+m)/(2*m)
end
imagesc(log10(draw+1))
%save 4000x4000t.30.rho.10.2.mat draw R X Y Z
print -dpng zoomindt22.png