Code covered by the BSD License

# MATLAB in Physics - Symbolic Computation and Differential Equations

### Matt McDonnell (view profile)

19 Feb 2009 (Updated )

The fourth lecture in a series on using MATLAB in undergraduate physics courses.

LorentzAttractor
```classdef LorentzAttractor < handle
%LORENTZATTRACTOR : simulation of the Lorentz Attractor chaotic system
%   The Lorentz attractor describes a model of the atmosphere that can
%   become chaotic in certain parameter regimes.  This class
%   encapsulates the simulation parameters and has methods to run the
%   simulation and plot the results.

%   Copyright 2008-2009 The MathWorks, Inc.
%   \$Revision: 35 \$  \$Date: 2009-05-29 15:27:34 +0100 (Fri, 29 May 2009) \$

properties
Sigma
Beta
Rho
MaxTime
SampleCount
StateValues
TimeValues
X0
end

methods
function obj = LorentzAttractor(varargin)
% Parse the input arguments using an input parser object
p = inputParser;
% sigma, beta and rho are the model parameters
% Maximum simulation time
% Number of points to return from the simulation result
% Starting state
% Parse the input arguments
p.parse(varargin{:});
% Store the parameters as properties of the object
obj.Sigma = p.Results.sigma;
obj.Beta = p.Results.beta;
obj.Rho = p.Results.rho;
obj.MaxTime = p.Results.tmax;
obj.SampleCount = p.Results.npts;
obj.X0 = p.Results.x0;
end

function simulate(obj)
% Perform the simulation
[obj.TimeValues,obj.StateValues] = ode45(...
@obj.propagate,[0 obj.MaxTime],obj.X0);
end

function y = propagate(obj,t,x) %#ok<INUSL>
% Equations of motion.  Time argument is unused for this system
% as the equations of motion are time independent, however the
% ODE solver requires equations of motion of the form
% f = f(t,x)
y = zeros(size(x));
y(1) = obj.Sigma.*(x(2) - x(1));
y(2) = x(1).*(obj.Rho - x(3)) - x(2);
y(3) = x(1).*x(2) - obj.Beta .* x(3);
end

function plot3(obj,varargin)
% Plot the system evolution in three dimensions.  If the
% simulation has not been performed yet then perform it here.
if isempty(obj.StateValues)
obj.simulate();
end
% Plot the evolution path
plot3(obj.StateValues(:,1),...
obj.StateValues(:,2),...
obj.StateValues(:,3),varargin{:});
% Show the starting point
startingPoint = line(obj.X0(1),...
obj.X0(2),...
obj.X0(3));
set(startingPoint,'LineStyle','none','Marker','o',...
'MarkerSize',8,'MarkerFaceColor','k');
% Set a nice view angle
view(45, 45);
end

function plot2D(obj,columnIndexes,varargin)
% Plot a 2D slice of the evolution, simulating the motion if
% the simulation has not been run already.
% columnIndexes is a 2 element vector of the dimensions to plot
% ie x=1, y=2, z=3
if isempty(obj.StateValues)
obj.simulate();
end

% Default to plotting y vs x if columnIndexes is unspecified
if nargin<2
columnIndexes = [1 2];
end

% Do the plot
plot(obj.X0(columnIndexes(1)),...
obj.X0(columnIndexes(2)),...
'k.',...
obj.StateValues(:,columnIndexes(1)),...
obj.StateValues(:,columnIndexes(2)),...
varargin{:});
end

function plotXY(obj,varargin)
% Plot the evolution in the XY plane
obj.plot2D([1,2],varargin{:});
end

function plotYZ(obj,varargin)
% Plot the evolution in the YZ plane
obj.plot2D([2,3],varargin{:});
end

function plotZX(obj,varargin)
% Plot the evolution in the ZX plane
obj.plot2D([3,1],varargin{:});
end
end

end

```