Differential Equations with Maple

The "Lorenz equation"

restart: with(DEtools):

An important problem in meteorology and other applications of fluid dynamics concerns

the motion of a layer of fluid, such as the earth's atmosphere, that is warmer at the bottom

than at the top.

If the vertical temperature difference NiMlJkRlbHRhRw==T is small then there is a linear variation of

temperature with altitude, but no significant motion of the fluid layer itself.

However if NiMlJkRlbHRhRw==T is large enough the warmer air begins to rise, displacing the cooler air

above it and a steady convective motion results.

If NiMlJkRlbHRhRw==T increases still further the steady convective flow breaks up and a more complex

turbulent motion is the result.

Edward Lorenz described this phenomenon around 1960 by a system of three first

order ODEs :

NiMvLSUlZGlmZkc2JC0lInhHNiMlInRHRiosJiomJSJhRyIiIkYnRi4hIiIqJkYtRi4tJSJ5R0YpRi5GLg==

NiMvLSUlZGlmZkc2JC0lInlHNiMlInRHRiosKComLSUiekdGKSIiIi0lInhHRilGLyEiIiomJSJyR0YvRjBGL0YvRidGMg==

NiMvLSUlZGlmZkc2JC0lInpHNiMlInRHRiosJiomLSUieUdGKSIiIi0lInhHRilGL0YvKiYlImJHRi9GJ0YvISIi

where a,b and r are constants.

In this sheet we'll investigate this system a little further.

dvx := diff(x(t),t) = -a * x(t) + a * y(t);

dvy := diff(y(t),t) = -z(t) * x(t) + r * x(t) - y(t);

dvz := diff(z(t),t) = y(t) * x(t) - b * z(t);

Let's assume the following values for the parameters a, b and r :

a:=10; b:=8/3; r:=28.0;

These values produce the wellknown picture for the Lorenz attractor.

deqns:=[dvx,dvy,dvz]:

vars:={x(t),y(t),z(t)}:

inits:=[[x(0)=5,y(0)=5,z(0)=5]]:

DEplot3d(deqns,vars, t=0..100, inits, stepsize=.02, x=-20..20, y=-25..25, z=0..50);

As you can see this output is not very clear : let's use some Maple tricks to rotate the result

somewhat and give the output a little more color :

DEplot3d(deqns,vars, t=0..100, inits, stepsize=.02, x=-20..20, y=-25..25, z=0..50, linecolour=sin(t*Pi/3), thickness = 1, orientation = [-40,80]);

If you want to investigate the projection of the trajectory on a certain plane (e.g. the

(x,y)- plane) you may use the following command :

phaseportrait(deqns,vars,t=0..20,inits,stepsize=0.02, scene=[x(t),y(t)],linecolour=blue);

Exercise 1 )

Now try to visualize the [x(t),z(t))] and [y(t),z(t)] planes by yourself.

To view e.g. x(t) as a function of t you can use the following command :

phaseportrait(deqns,vars,t=0..20,inits,stepsize=0.02, scene=[t,x(t)],linecolour=green);

Exercise 2 )

The phenomenon of chaotic behaviour only occurs for values of r being large enough.

Investigate the behaviour of x(t) as function of t for values of r being respectively :

r = 1

r = 10

r = 25

and tell whether the behaviour is chaotic or not. Use the theory in your textbook

(e.g. page 532 - 538 in Boyce & Diprima, 7th edition) to check your answers.

Note: you have to redefine the equations after changing the value of r to ensure

the new value of r is being used.

r := ??: deqns:=[dvx,dvy,dvz]:

phaseportrait(deqns,vars,t=0..20,inits,stepsize=0.02,scene=[t,x(t)],linecolour=green);

Exercise 3 )

Now restore the original parameters for the Lorenz equations from above and try to

start from 3 different initial points that you may choose for yourself.

Describe the effect of the initial point on the "steady-state" solution.

r := 28: deqns:=[dvx,dvy,dvz]:

inits:=[[x(0)=??,y(0)=??,z(0)=??]]:

DEplot3d(deqns,vars, t=0..100, inits, stepsize=.02, x=-20..20, y=-25..25, z=0..50, linecolour=sin(t*Pi/3), thickness = 1, orientation = [-40,80]);

version 1.3, KL, Nov 2006