Code covered by the BSD License  

Highlights from
Solution of Differential Equations with MATLAB & Simulink: Lorenz Attractor Case Study

image thumbnail

Solution of Differential Equations with MATLAB & Simulink: Lorenz Attractor Case Study



11 Aug 2010 (Updated )

Simulink design pattern for solving differential equations, visualize results in MATLAB graphics

lorenz_graphs_level_1(t,x,u,flag, x_initial)
function [sys,x0,str,ts,simStateCompliance] = lorenz_graphs_level_1(t,x,u,flag, x_initial)
% Level 1 S-Functions are no longer supported. Use this template below to
% understand any legacy code implemented in this template. 
% The goal of this code is to make the S-function behave like a Simulink
% block. During simulation of a model, the Simulink engine invokes this function with different flag values so that
% sub functions associated with those values get called. The S-function
% performs the task and returns the results in an output vector. For
% example, flag 0 will invoke the initailization of this block during the
% initialization of the entire model
% Function parameter definition:
% Inputs
%   t	Current time
%   x	State vector
%   u	Input vector
%  flag	Integer value that indicates the task to be performed by the S-function

% Outputs
% sys

% Depending on the value of the flag variable, the following functions will
% be called during the execution.
%	Flag value      Method               What it does (in brief)
%        0          mdlInitializeSizes	 Initialize sys
%        1	        mdlDerivatives	     Calculates the derivatives of the continuous state variables.
%        2	        mdlUpdate	         Updates discrete states, sample times, and major time step requirements.
%        3	        mdlOutputs	         Calculates the outputs of the S-function.
%        4	        mdlOutputs           Updates the run-time object NextTimeHit property	
%        9          mdlTerminate         Perform any necessary end-of-simulation tasks
% Additional block parameters are appended at the end of the function
% prototype declaration. Here we want the initial states of the
% differential equation that is in sync with the initial conditions used in
% the integrator block. 
% Remember that if you have extra parameters that are being input to the
% block, then the function prototype will have the extra variables in the
% same order appended at the end. i.e. x_initial. For passing these
% parameters around, you need to pass them as arguments to the rest of the
% function. This is not the case for Level-2 S-functions that do not
% require such explicit modification.

switch flag,
  case 0,
     % Modify the InitialSizes Function so that the initial conditions can
     % be passed in the call
  case 2,
  case { 1, 3, 4, 9 },
    sys = [];

function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes(x_initial)

% Declare an uninitialized simsez structure that will contain information
% about the system
sizes = simsizes;

% No. of continuous states
sizes.NumContStates  = 0;

% No. of discrete states
sizes.NumDiscStates  = 3;

% No. of outputs
sizes.NumOutputs     = 0;

% No. of inputs
sizes.NumInputs      = 3;

% Flag for direct feedthrough
sizes.DirFeedthrough = 0;

% No. of sample times
sizes.NumSampleTimes =1;

% Pass information from sizes structure to the sys structure which the
% Simulink engine will use
sys = simsizes(sizes);

% Initialize the initial conditions
x0  = [x_initial]';

% Reserved for future use, initialize it to []
str = [];

% For the sample time to inherit from the model use the following. 
% If you want it to run every 0.25 seconds (discrete sample time) starting at 0.1 seconds after the simulation start time, set ts to [0.25 0.1].
% This will slow the simulation down but will improve the quality of the
% plot. This is a tradeoff in the simulation

% If you want to inherit from the model itself
%ts =[-1 0];

% Visualization in a MATLAB GUI
% Initialization code
set(q, 'Name', 'Level-2 Lorenz Oscillator Visualization', ...
    'Color', [1 1 1]);
set(gca,'Units', 'normalized');
title('Position (x,y,z)');
hold on;
grid on
axis([-50 50 -100 100 0 200])
axis vis3d % freezes the aspect ratio to prevent auto strech to fill which can introduce jarring animation
set(gca, 'xtick',[-50:25:50], 'ytick',[-100:50:100], 'ztick',[0:50:200]);
set(gca,'xColor', [1 0 1], 'ycolor', [0 0 0], 'zcolor',[0 0 1])

simStateCompliance = 'DefaultSimState';

function sys=mdlUpdate(t,x,u) %#ok<INUSL>
% Connect the current computed point with the past connected point
line([u(1,:) x(1, :)], [u(2,:) x(2,:)], [u(3,:) x(3,:)],'LineWidth',2,...
    'MarkerFaceColor',[.49 1 .63],...
    'MarkerSize',12, 'EraseMode', 'none');
% Plot an asterisk
a=plot3(u(1,:), u(2,:), u(3,:),'r.', 'EraseMode', 'none', 'MarkerSize', 22);

% Comet command will add a comet to the tip
%comet3(u(1,:), u(2,:), u(3, :), 0.01);

% Rotate the camera
camorbit(0.5,0.5, 'data', [1 1 1])

% Update the points plotting

% The state is the output which is stored for the next time step. This is
% for the internal state since we are calling mdlUpdate

Contact us