function lorenz

% the above is required, if function definition is incldued within the
% program. Do not remove the above statement
%
% an example based on 
% R. C. Hilborn, "Chaos and Nonlinear Dynamics - an introduction for
% scientists and engineers", second edition, Oxford University Press, 2000.
% page 27-37.

% ========================================================================
%
% The matlab code for Chaos Theory - 
%                 Nonlinear Statistical Modeling of Speech
% date : July 15, 2005 Friday.
% IES: http://www.cavs.msstate.edu/hse/ies/projects/nsf_nonlinear
%
% ========================================================================

% setting up the environment
% clearing all the variables, clearing the command window, closing all the
% figures (plots)
%
close all;
clear all;
clc;

% In this example we change the error tolerances using the odeset command 
% and solve on a time interval [0 12] with 
% an initial condition vector [0 1 0] at time 0.
%
xinitial = 0;
yinitial = 1;
zinitial = 0;

tstart = 0;
tend = 12;

options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
[T,Y] = ode45(@rigid,[tstart tend],[xinitial yinitial zinitial], options);

% solve the problem using ODE45
% plotting the results for X, Y, Z versus time,
% and state-space projection. 
% All the graphs are put on a single graph page
%
figure;
subplot(221); plot(T,Y(:,1),'-'); 
xlim([tstart tend]); ylim([min(Y(:,1)) max(Y(:,1)) + 0.1]);
title('Lorenz eqn - X vs. time');

subplot(222); plot(T,Y(:,2),'-.');
xlim([tstart tend]); ylim([min(Y(:,2)) max(Y(:,2))]);
title('Lorenz eqn - Y vs. time');

subplot(223); plot(T,Y(:,3),'.');
xlim([tstart tend]); ylim([min(Y(:,3)) max(Y(:,3)) + 0.1]);
title('Lorenz eqn - Z vs. time');

subplot(224); plot(Y(:,1),Y(:,3),'.');
xlim([min(Y(:,1)) max(Y(:,1)) + 0.1 ]);
ylim([min(Y(:,3)) max(Y(:,3)) + 0.1]);
title('Lorenz eqn - Z vs. X - state space Projection');

figure;
plot3(Y(:,1), Y(:,2), Y(:,3));
grid on;
title('Phase state representation - Lorenz Equations');
% --------------------------------------------------------------------------
function dy = rigid(t,y)

% the lorenz model equations are given below
% the three equations are:
% dx / dt = p(Y - X)          ...........eqn(1)
% dy / dt = -XZ + rX - Y      ...........eqn(2)
% dz / dt = -XY -bZ           ...........eqn(3)
% p = Prandtl number, r = Rayleigh number, b = vertical to horizontal ratio
%
r = 0.5;
p = 10;
b = 8.0/3;

dy = zeros(3,1);    
dy(1) = -1 * p * y(1) + p * y(2); 
dy(2) = (-r + y(3)) * y(1) -1 * y(2);
dy(3) = y(2) - b * y(3);


