Trisha’s Blog!

Block 3: The Lorenz Attractor

In this block we are working on systems of Non-Linear Differential Equations, The Edward Lorenz’s equations, and the Lorenz attractor. When dealing with biological and physical phenomena it can be hard to find solution curves to your differential equation. In fact there are very few exact answers, so instead we have to find many solution curves to get and approximation of the answer. We have already practiced these methods with the Euler’s and Runge-Kutta method. The problem is Euler creates a very large margin or error due to the step size. Runge-Kutta, however, creates a much more accurate solution curve and can help us when working with the Lorenz method. Matlab has a built in program that demonstrates the Lorenz attractor and how it works. Edward Lorenz created a simplified version of the convection rolls that arise in the atmosphere. He simplified the equation into 3 separate equations:

\frac{dx}{dt} = \sigma (y - x)

\frac{dy}{dt} = x (\rho - z) - y

\frac{dz}{dt} = xy - \beta z

All of these variables; \sigma,\rho, \beta are positive. The purpose of this activity is to show how Edward Lorenz rediscovered the chaos theory. He showed how the sensitive the equations are to initial conditions, known as the butterfly effect. \sigma  and \rho  are related to physical properties of a fluid in motion. With problems like this it is most common that \sigma  and \beta  have values of \sigma = 10, \beta = 8/3.  Changing the value of \rho is what gives you the different solution curves to these systems of equations.

Now I will demonstrate the differences between how Euler solves a system of D.E’s versus how Lorenz does. With Euler we can show as a single graph like we did in block one or we can tell Matlab that y is a column vector with as many components (= rows) as there are unknown functions. This will give us a 3-D graph of our answer. I will demonstrate both.

To create a graph using the Euler method we first need to write out a code for matlab. A basic Euler code saved as an m-file will be fine, I used the one given on the class website, it is saved as Euler.m:

% Euler’s method for a single differential equation.

% The differential equation is dy/dt = f(t,y).

function [ t, y ] = euler ( f, t_range, y_initial, nstep )

% We set a range of time values, from t(1) to t(2).

t(1) = t_range(1);

% We define dt by dividing the time range into an equal number of specified steps.

dt = ( t_range(2) – t_range(1) ) / nstep;

% We set the initial value of y at the beginning time t(1).

y(1) = y_initial;

% We use Euler’s method to update the value of y at new time steps.

%  “feval” is used instead of “eval” because we are passing the name of f to the program.

for i = 1 : nstep
t(i+1) = t(i) + dt;
y(i+1) = y(i) + dt * feval ( f, t(i), y(i) );
end                                                                                                        plot(t,y)                               

Next we need a second m-file that describes this function. Again I will use the one given on the class website it is saved as Test_example.m:                                                                                  

function yprime = test_example ( t, y )
yprime = y*(2./t-1);

Now we enter on Matlab:

>> y_init = 0.1;
>> [ t, y ] = euler ( ‘test_example’, [ 0.1, 9.0 ], y_init, 200 );

This code then produces a graph:

euler-lorenz

Now that we have this graph we can modify the Euler m-file so that it can calculate a numerical approximation to the lorenz equations. To do this we basically use the same code as above but change all of the “y(1) or y(i+1)” to “y(:,1) and y(:,i+1)” and at the bottom we add: plot3(y(1, :),y(2, :),y(3, :)). Adding the “:,” to the y values just tells the computer that  y is a column vector with as many components as there are unknown functions. I will save this as euler_system.m.

Next we create a another m-file that demonstrates the lorenz function. I am using the one given on the class website. This code is saved as Lorenz_system.m.

function yprime = lorenz_system ( t, y )
yprime = [ 10.0* (y(2)-y(1)); y(1)*(28.0-y(3))-y(2);y(1)*y(2)-8*y(3)/3 ];

Finally, we go to Matlab and enter in the codes listed below to run the Euler method of the Lorenz function. These codes then produce the graph below.

>> y_init = [ rand(); rand()rand() ];

>> [ t, y ] = euler_system ( ‘lorenz_system’, [ 0.0, 20.0 ], y_init, 1000 );euler-lorenz2

Now we can use Matlab’s built-in program Ode45 to create possible solution curves to the Lorenz system. This process will require two more m-file’s which I have listed below. These were given on the class website. The first is saved as g.m and the second as lorenz_demo.m 

g.m:

function xdot = g(t,x)

xdot = zeros(3,1);

sig = 10.0;

rho = 28.0;

bet = 8.0/3.0;

xdot(1) = sig*(x(2)-x(1));

xdot(2) = rho*x(1)-x(2)-x(1)*x(3);

xdot(3) = x(1)*x(2)-bet*x(3);

 

Lorenz_demo:

function lorenz_demo(time)

% Usage: lorenz_demo(time)

% time=end point of time interval

% This function integrates the lorenz attractor

% from t=0 to t=time

[t,x] = ode45(’g’,[0 time],[1;2;3]);

disp(’press any key to continue …’)

pause

plot3(x(:,1),x(:,2),x(:,3))

print -deps lorenz.eps

Finally, with these two new m-files we can go back into Matlab and enter the codes: >>Lorenz_demo(200). Which then produces the graph shown below.lorenzdemo

Rossler system:

Now that we have explored the Lorenz Equations and can see that they are in fact chaotic, we can explore new ways to show this. I think The Rossler system would be a good start. When demonstrating this system there are set values already established for this system.

a=.1

b=.1

c=14

Now we can set-up an m-file:

function yprime=rossler_system1(t,y)
yprime = [-y(2)-y(3);y(1)+0.1*y(2);0.1+(y(1)-14.0)];

In Matlab we enter

>> y_init = [ rand(); rand(); rand() ];
>> [ t, y ] = euler_system ( ‘rossler_system1′, [ 0.0, 20.0 ], y_init, 1000 );

which gives us the graph

b3-rossler1-3d

Now we can use the Rung-Kutta method to show this graph in more detail. First, though we need to create two more m-files:

-function xdot=r(t,x)
xdot=zeros(3,1);
alp=0.1;
bet=0.1;
gam=14.0;
xdot(1)=-x(2)-x(3);
xdot(2)=x(1)+ alp*x(2);
xdot(3)=bet+x(3)*(x(1)-gam);

-function rossler_demo(time)
[t,x] = ode45(’r’,[0 time],[1;2;3]);
plot3(x(:,1),x(:,2),x(:,3)

Then in Matlab we enter 

>>rossler_demo(200)

And we get this pretty graph:

b3-rossler1-demo

 

Now we can see how the Rossler system changes when we make the slightest adjustments to its values. I changed the values of the initial rossler system to

a=

b=

c=

Using the same m-file as before just adjusting the values to mach our new ones.

function yprime=rossler_system1(t,y)
yprime = [-y(2)-y(3);y(1)+0.1*y(2);0.1+(y(1)-14.0)];

Next we execute the same command in Matlab

>> y_init = [ rand(); rand(); rand() ];
>> [ t, y ] = euler_system ( ‘rossler_system1′, [ 0.0, 20.0 ], y_init, 1000 );

and we get this graph below.

Rossler EX1

As you can see the system is indeed chaotic just as the Lorenz system. Unfortunately I am experiencing difficulty with Matlab and demonstrating this graph using the Rung-Kutta method. I have checked over my m-files many times and can still see no reason why Matlab is not producing the graph for me.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

%d bloggers like this: