Convection

currents
Convection currents on earth mantle produce tectonic displacement
 Foundation of chaos theory
Lorenz did not has a title for a lecture to be
given in 1972 at the 139th congress of the
The butterfly effect
American Association for the Advancement of Foundation of the
Science. Philip Merilees, the chairperson, Popular Fascination
presented it as: “Does the flap of a butterfly’s
wings in Brazil set off a tornado in Texas?”.
produced by chaos
Theory
 Lorenz Equations
dx
 ay  ax
dt
dy
 bx  y  zx
dt
dz
 xy  cz
dt

t : time

¿ ?
x(t) : position
y(t) : temperature
z(t) : “distortion of T”
a, b, c : Parameters
Computer solutions to differential equations

x  f ( x), x(0)  x0
We ride on the phase point, flowing according to the
vector field.

If we begin at t=0 in x(0)=x0 with a velocity f(x0), after a

very small period of time t we would advance a
distance f(x0)t so now we will be at x(t)=x1

x1 = x0 + f(x0)t
Once again
x2 = x1 + f(x1)t
And, after several iterations,
xn+1 = xn + f(xn)t
Euler Method
 d
x(t )  sin  x(t )  , x(0)  x0 through Euler
dt
6

x0 = pi*[-2 -1.9 -3/2 -1 -1/2 -0.1 0 0.1

1/4 1/2 1 3/2 1.9 2];
for i=1:14
x(1) = x0(i); 4
for j=2:10000
x(j) = x(j-1) + 0.0004*sin(x(j-1));
end
plot(0.0004*(0:9999),x)
hold on
end 2

x2
0

x1 -2
-2 - 0  2

-4

-6
 d
x(t )  sin  x(t )  , x(0)  x0 through Euler
dt
• Just estimate the
derivative at the left
x(1) = 0.01; y(1) = 0.01;
for j=2:20000 end of the interval
if j==200*floor(j/200)
y(1+j/200) = y(j/200) + 0.08*sin(y(j/200));
(tn,tn+1)
end
x(j) = x(j-1) + 0.0004*sin(x(j-1));
• Better to use the
end average of the
plot(0.0004*(0:19999),x,'b-',0.08*(0:100),y,'r--')
derivative within the
4
interval
20000 pasos d  xn  f ( xn )t
100 pasos xn1  xn  12  f ( xn )  f (d )  t
3
estado, x(t)

0
0 2 4 6 8
tiempo, t
d
x(t )  sin  x(t )  , x(0)  x0 Through 2nd order Euler
dt
x(1) = 0.01; y(1) = 0.01;
for j=2:20000
if j==200*floor(j/200)
d = y(j/200) + 0.08*sin(y(j/200));
y(1+j/200) = y(j/200) + 0.04*(sin(y(j/200)) + sin(d));
end
x(j) = x(j-1) + 0.0004*sin(x(j-1));
end
plot(0.0004*(0:19999),x,'b-',0.08*(0:100),y,'r--')
legend('20000 pasos','100 pasos')
xlabel('tiempo, t')
ylabel('estado, x(t)')
First order Euler
3.5 Second order Euler
20000 pasos
3
100 pasos
2.5
estado, x(t)

1.5

0.5

0
0 1 2 3 4 5 6 7 8
tiempo, t
 4th order Runge-Kutta method
k1  f ( xn )t
k2  f  xn  k1 / 2  t
k3  f  xn  k2 / 2  t
k4  f  xn  k3  t
xn 1  xn  16  k1  2k2  2k3  k4 
4
20000 pasos
20 pasos
3
estado, x(t)

0
0 2 4 6 8
tiempo, t
 Multidimensional 4th order Runge-Kutta method
k1  f ( xn )t

k2  f xn  k1 / 2 t 
k3  f x n k / 2  t
2

k4  f x n k  t
3

xn 1  xn   k  2k  2k
1
6 1 2 3  k4 
function x = AtractorDeLorenz
x(:,1) = [1 2 3]';
for j=2:10000
dx k1 = 0.01*ecuacion(x(:,j-1));
 ay  ax k2 = 0.01*ecuacion(x(:,j-1) + k1/2);
dt k3 = 0.01*ecuacion(x(:,j-1) + k2/2);
k4 = 0.01*ecuacion(x(:,j-1) + k3);
x(:,j) = x(:,j-1) + (k1 + 2*k2 + 2*k3 + k4)/6;
dy
 bx  y  zx end
x = x';
dt plot3(x(:,1),x(:,2),x(:,3));

dz function x = ecuacion(x)
 xy  cz x = [10*(x(2)-x(1)); 28*x(1)-x(2)-x(1)*x(3); x(1)*x(2)-(8/3)*x(3)];

dt
Of course, in matlab Runge-Kutta is already implemented!
g = inline(‘[10*(x(2)-x(1)); 28*x(1)-x(2)-x(1)*x(3); x(1)*x(2)-(8/3)*x(3)]‘,'t','x');
[t,x] = ode45(g,[0:0.01:100],[1;2;3]);
plot(t,x);

Time evolution: ¡Chaos!

20

10

x(t) 0

-10

-20
0 10 20 30 40 50 60 70 80 90 100

30

20

10

y(t) 0

-10

-20

-30
0 10 20 30 40 50 60 70 80 90 100

50

40

30
z(t) 20

10

0
0 10 20 30 40 50 60 70 80 90 100

Something more interesting?

 g = inline(‘[10*(x(2)-x(1)); 28*x(1)-x(2)-x(1)*x(3); x(1)*x(2)-(8/3)*x(3)]','t','x');
[t,x] = ode45(g,[0:0.01:100],[1;2;3]);
plot3(x(:,1),x(:,2),x(:,3));

50

45
Phase Space
40

35

30

25

20

15

10

0
-30
-20
-10
0
10
20 10 15 20
30 -5 0 5
-20 -15 -10
 Homework #6, part III

• Consider the following 3-dimensional differential equation (Rossler equation):

d
x  ( y  z )
dt
d
y  x  0.2 y
dt
d
z  0.2  z ( x  5.7)
dt
• Use the Runge-Kutta numerical integration method to plot x(t) vs t, y(t) vs t and z(t) vs t for two slightly
different initial conditions.
• Use the Runge-Kutta numerical integration method to draw a 3-dimensional plot of the trajectory (a
strange attractor)
• Discuss the results
 A (very) brief introduction to chaos control
OGY : Ott E., Grebogi C., Yorke J. "Controlling Chaos" Physics Review Letters, 1990

xn   xn1 1  xn1  , 0  xn  1
• In its natural state,  corresponds to a chaotic regime, e.g.  = 3.78
• For some reason, we cannot move  to an stable regime
• We can only produce very small perturbations at particular instants
of time, where the system remains in a chaotic regime
• We take advantage of the existence of a fixed point (unstable, of
course) at x* = ( - 1)/ 
• When x(t) gets close to x*, we adjust 

x(1) = 0.5; % Punto inicial

ro = 3.78; % Parámetro en régimen caótico
r = ro; % Nuestro parámetro
for i=2:2500
x(i) = r*x(i-1)*(1-x(i-1)); % Evalúa el siguiente estado
rp = 1/(1-x(i)); % Y, si está muy cerca del punto fijo
if abs(rp-ro)<0.01 % hace un ligerísimo ajuste al parámetro
r = rp;
end
end
plot(x)
 1

x(t) 0.5

0
0 500 1000 1500 2000 2500

3.8

(t) 3.79

3.78
0 500 1000 1500 2000 2500

In the computer, the system will remain at

the equilibrium point but, in reality, noise or 1

Imprecision will take it out of control. 0.8

0.6

0.4

x(1) = 0.5; % Punto inicial 0.2

ro = 3.78; % Parámetro en régimen caótico
0
r(1) = ro; % Nuestro parámetro 0 500 1000 1500 2000 2500
for i=2:2500
x(i) = r(i-1)*x(i-1)*(1-x(i-1)) + 0.0001*randn(1); % Evalúa el siguiente estado 3.795
rp = 1/(1-x(i)); % Y, si está muy cerca del punto fijo
if abs(rp-ro)<0.01 % hace un ligerísimo ajuste al parámetro 3.79
r(i) = rp;
3.785
else
r(i) = r(i-1); 3.78
end
end 3.775
0 500 1000 1500 2000 2500
subplot(211); plot(x)
subplot(212); plot(r,'r')
 1

0.8

0.6

0.4

0.2

0
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000

3.795

3.79

3.785

3.78

3.775

3.77
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000

If noise is too high, the system can get out of control every now and then
 Period-2 periodic orbit
1

0.8

0.6

0.4

0.2

0
0 50 100 150

3.784

3.783

3.782

3.781

3.78

3.779
0 50 100 150

x(1) = 0.5;
ro = 3.78;
r(1) = ro;
x(2) = r(1)*x(1)*(1-x(1));
r(2) = ro • Chaos control in ventricular fibrillation
for i=3:5000 • Chaos control in brain activity
rp = x(i-2)/(x(i-1)*(1-x(i-1)));
if abs(rp-ro)<0.01 associated with epileptic attack
r(i) = rp;
else
r(i) = r(i-1);
end
x(i) = r(i)*x(i-1)*(1-x(i-1));
end
subplot(211); plot(x)
subplot(212); plot(r,'r')