pplane7.m

matlab of linear system

function output = pplane7(action,input1,input2,input3)

% pplane7  is an interactive tool for studying planar autonomous systems of
%	differential equations.  When pplane7 is executed, a pplane7 Setup
%	window is opened.  The user may enter the differential
%	equation and specify a display window using the interactive
%	controls in the Setup window.  Up to 4 parameters may also be 
%	specified.  In addition the user is give a choice of the type of 
%	field displayed and the number of field points.
%
%	When the Proceed button is pressed on the Setup window, the pplane7
%	Display window is opened, and a field of the type requested is 
%	displayed for the system.  When the mouse button is depressed in 
%	the pplane7 Display window, the solution to the system with that
%	initial condition is calculated and plotted.
%
%	Other options are available in the menus.  These are
%	fairly self explanatory.  
%
%	This is version 6.0, and will only run on version 6 of MATLAB.

%  Copywright (c)  1995, 1997, 1998, 1999, 2000, 2001, 2002, 2003 
%                  John C. Polking, Rice University
%                  Last Modified: April 24, 2003

startstr = 'pplane7';
if nargin < 1
  action ='initialize';
end



if get(0,'screendepth') == 1
   style = 'bw';
end

if strcmp(action,'initialize')
   
   % First we make sure that there is no other copy of pplane7
   % running, since this causes problems.	
   pph = findobj('tag','pplane7');
   if ~isempty(pph);
      qstring = {'There are some pplane7 figures open although they may be invisible.  ';...
            'What do you want to do?'};
      tstring = 'Only one copy of pplane7 can be open at one time.';
      b1str = 'Restart pplane7.';
      b2str = 'Just close those figures.';
      b3str = 'Do nothing.';
      answer = questdlg(qstring,tstring,b1str,b2str,b3str,b1str);
      if strcmp(answer,b1str)
         delete(pph);
         %pplane7;return
         eval(startstr);return
      elseif strcmp(answer,b2str)
         delete(pph);return
      else
         return
      end 
   end  % if ~isempty(pph);
   
   % Make sure tempdir is on the MATLABPATH.  We want to be sure that we
   % return the path to its current position when we exit.
   
   p = path;
   tmpdir = tempdir;
   ll = length(tmpdir);
   tmpdir = tmpdir(1:ll-1);
   ud.remtd = 0;
   if isempty(findstr(tmpdir,p))
     ud.remtd = 1;
     addpath(tempdir)
   end
   
   % Next we look for old files created by pplane7.
   % First in the current directory.
   
   oldfiles = dir('pptp*.m');
   ofiles = cell(0,1);
   kk = zeros(0,1);
   for k = 1:length(oldfiles)
      fn = oldfiles(k).name;
      fid = fopen(fn,'r');
      ll = fgetl(fid);
      ll = fgetl(fid);
      ll = fgetl(fid);
      fclose(fid);
      if strcmp(ll,'%% Created by pplane7')
         delete(fn)
      else
	l= length(ofiles);
	ofiles{l+1} = fn;
      end
   end
   
   %Then in the temp directory.
   
   oldfiles = dir([tempdir,'pptp*.m']);
   for k = 1:length(oldfiles)
      fn = [tempdir,oldfiles(k).name];
      fid = fopen(fn,'r');
      ll = fgetl(fid);
      ll = fgetl(fid);
      ll = fgetl(fid);
      fclose(fid);
      if strcmp(ll,'%% Created by pplane7')
         delete(fn)
      else
	l= length(ofiles);
	ofiles{l+1} = fn;
      end
   end
   lll = length(ofiles);
   if lll >0
     if lll == 1
       astr = 'The file';
       bstr = 'If so it can be safely deleted.';
     else
       astr = 'The files';
       bstr = 'If so they can be safely deleted.';
     end
     fprintf([astr,'\n']);
     for j = 1:lll
       fn = ofiles{j};
       disp(['     ',fn]);
     end
     fprintf('may have been created by an older version of PPLANE.\n');
     fprintf(bstr,'\n\n');
   end
   
   style = 'white';
   ppdir = pwd;
   ssize = get(0,'defaultaxesfontsize');
   npts = 20;
   solver = 'Dormand Prince';
   tol = 1e-4;
   stepsize = 0.1;
   if exist('ppstart','file')
     H = ppstart;
     if ~isempty(H)
       if isfield(H,'style')
	 style = H.style;
       end
       if isfield(H,'size')
	 ssize = H.size;
       end
       if isfield(H,'npts')
	 npts = H.npts;
       end
      if isfield(H,'solver')
	solver = H.solver;  
      end
      if isfield(H,'ppdir')
	ppdir = H.ppdir;  
      end
      if isfield(H,'stepsize')
	stepsize = H.stepsize;  
      end
      if isfield(H,'tolerance')
	tol = H.tolerance;  
      end
    end
  end
  if get(0,'screendepth') == 1
    style = 'bw';
  end
  
  ud.ssize = ssize;
  ud.ppdir = ppdir;
  comp = computer;
  if strcmp(comp,'PCWIN')
    ud.fontsize = 0.8*ud.ssize;
  else
    ud.fontsize = ud.ssize;
  end
   
  % Set up for the menu of systems.
  
  system.name = 'default system';
  system.xvar = 'x';
  system.yvar = 'y';
  system.xder = ' 2*x - y + 3*(x^2-y^2) + 2*x*y';
  system.yder = ' x - 3*y - 3*(x^2-y^2) + 3*x*y';
  system.pname = {};
  system.pval = {};
  system.fieldtype = 'arrows';
  system.npts = npts;
  system.wind = [-2 4 -4 2];
  
  system(2).name = 'linear system';
  system(2).xvar = 'x';
  system(2).yvar = 'y';
  system(2).xder = ' A*x + B*y';
  system(2).yder = ' C*x + D*y';
  system(2).pname = {'A', 'B', 'C', 'D', '', ''};
  system(2).pval = {'2', '2', '-2',  '-3'};
  system(2).fieldtype = 'arrows';
  system(2).npts = npts;
  system(2).wind = [-5 5 -5 5];
  
  system(3).name = 'vibrating spring';
  system(3).xvar = 'x';
  system(3).yvar = 'v';
  system(3).xder = ' v';
  system(3).yder = ' -(k*x + d*v)/m';
  system(3).pname = {'k', 'm', 'd', '', '', ''};
  system(3).pval = {'3', '1', '0'};
  system(3).fieldtype = 'arrows';
  system(3).npts = npts;
  system(3).wind = [-5 5 -5 5];
  
  system(4).name = 'pendulum';
  system(4).xvar = '\theta';
  system(4).yvar = '\omega';
  system(4).xder = ' \omega';
  system(4).yder = ' -sin(\theta) - D*\omega';
  system(4).pname = {'D', '', '', '', '', ''};
  system(4).pval = {'0'};
  system(4).fieldtype = 'arrows';
  system(4).npts = npts;
  system(4).wind = [-10 10 -4 4];
  
  system(5).name = 'predator prey';
  system(5).xvar = 'prey';
  system(5).yvar = 'predator';
  system(5).xder = ' (A - B*predator)*prey';
  system(5).yder = ' (D*prey - C)*predator';
  system(5).pname = {'A', 'B', 'C', 'D', '', ''};
  system(5).pval = {'0.4', '0.01', '0.3', '0.005'};
  system(5).fieldtype = 'arrows';
  system(5).npts = npts;
  system(5).wind = [0 120 0 80];
  
  system(6).name = 'competing species';
  system(6).xvar = 'x';
  system(6).yvar = 'y';
  system(6).xder = ' r*(1 - x - A*y)*x';
  system(6).yder = ' s*(1 - y - B*x)*y';
  system(6).pname = {'r', 's', 'A', 'B', '', ''};
  system(6).pval = {'0.4', '0.6', '5', '4'};
  system(6).fieldtype = 'nullclines';
  system(6).npts = npts;
  system(6).wind = [0 1 0 1];
   
  system(7).name = 'cooperative species';
  system(7).xvar = 'x';
  system(7).yvar = 'y';
  system(7).xder = ' r*(1 - x + A*y)*x';
  system(7).yder = ' s*(1 - y + B*x)*y';
  system(7).pname = {'r', 's', 'A', 'B', '', ''};
  system(7).pval = {'0.4', '0.6', '2', '0.3'};
  system(7).fieldtype = 'nullclines';
  system(7).npts = npts;
  system(7).wind = [0 9 0 5];
  
  system(8).name = 'van der Pol''s equation';
  system(8).xvar = 'x';
  system(8).yvar = 'y';
  system(8).xder = ' M*x - y - x^3';
  system(8).yder = ' x';
  system(8).pname = {'M', '', '', '', '', ''};
  system(8).pval = {'2'};
  system(8).fieldtype = 'arrows';
  system(8).npts = npts;
  system(8).wind = [-5 5 -5 5];
  
  system(9).name = 'Duffing''s equation';
  system(9).xvar = 'x';
  system(9).yvar = 'y';
  system(9).xder = ' y';
  system(9).yder = ' -(k*x + c*y + l*x^3)/m';
  system(9).pname = {'k', 'c', 'm', 'l', '', ''};
  system(9).pval = {'-1', '0', '1', '1'};
  system(9).fieldtype = 'arrows';
  system(9).npts = npts;
  system(9).wind = [-3 3 -3 3];
  
  system(10).name = 'square limit set';
  system(10).xvar = 'x';
  system(10).yvar = 'y';
  system(10).xder = ' (y + x/5)*(1-x^2)';
  system(10).yder = ' -x*(1-y^2)';
  system(10).pname = {'', '', '', '', '', ''};
  system(10).pval = {'', '', '', ''};
  system(10).fieldtype = 'arrows';
  system(10).npts = npts;
  system(10).wind = [-1.5 1.5 -1 1];
  
  
  ud.c = system(1);	% Changed values.
  pname = ud.c.pname;
  for kk = length(pname)+1:6
    pname{kk} = '';
  end
  pval = ud.c.pval;
  for kk = length(pval)+1:6
    pval{kk} = '';
  end
  ud.c.pname = pname;
  ud.c.pval = pval;
  ud.o = ud.c;	% Original values.
  % ud.h = system(1);	% This will be the handles in the
			% setup window.
			
  ud.style = style;
  ud.npts = npts;  
  ud.settings.magn = 1.25;
  ud.settings.refine = 8;
  ud.settings.tol = tol;
  ud.settings.solver = solver;
  ud.settings.stepsize = stepsize;
  ud.settings.speed = 100;
  ud.system = system;
  ud.solvers = {'ppdp45';
		'rk4';
		'ode45';
		'ode23';
		'ode113';
		'ode15s';
		'ode23s';
		'ode23t';
		'ode23tb'};
  
  switch style
   case 'black'
    color.temp = [1 0 0]; % red for temporary orbits
    color.orb = [1 1 0];  % yellow for orbits
    color.eqpt = [1 0 0];  % red for eq. pts.
			   %    color.arrows = [0 1 1]; % cyan for arrows
			   % color.arrows = [.5 .5 .9];  % purple for arrows
    color.arrows = .5*[1 1 1];  % gray for arrows
    color.narrows = .7*[1 1 1];  % gray for nullcline arrows
    color.tx = [1 1 0]; % yellow for xt plots & 3D plots
    color.ty = [1 0 0]; % red for yt plots
    color.ta = [1 0 0]; % red for axis plots
    color.sep = [.2,1,0]; % green for separatrices
    color.xcline = [1 1 1]; % white for xclines
    color.ycline = [1 0 1]; % magenta for yclines
    color.level = [1,.5,.5];
    
   case 'white'
    color.temp = [1 0 0]; % red for temporary orbits
    color.orb = [0 0 1];  % blue for orbits
    color.eqpt = [1 0 0];  % red for eq. pts.
    color.arrows = 0.7*[1 1 1]; % gray for arrows
    color.narrows = .4*[1 1 1];  % gray for nullcline arrows
    color.tx = [0 0 1]; % blue for xt plots & 3D plots
    color.ty = [1 0 0]; % red for yt plots
    color.ta = [1 0 0]; % red for axis plots
    color.sep = [0,1,0];% green for separatrices
    color.xcline = [1 0 .75]; % magenta for xclines
    color.ycline = [1 .5 0]; %  orange for yclines
    color.level = 0.8*[.9,.5,.8];
    
   case 'test'
    color.temp = [1 0 0]; % red for temporary orbits
    color.orb = [0 0 1];  % blue for orbits
    color.arrows = .7*[1 1 1]; % gray for arrows
    color.eqpt = [1 0 0];  % red for eq. pts.
    color.narrows = .4*[1 1 1];  % gray for nullcline arrows
    color.tx = [0 0 1]; % blue for xt plots & 3D plots
    color.ty = [1 0 0]; % red for yt plots
    color.ta = [1 0 0]; % red for axis plots
    color.sep = [0,1,0];% green for separatrices
    color.xcline = [1 0 .75]; % magenta for xclines
    color.ycline = [1 .5 0]; %  orange for yclines
    color.level = [1,.5,.5];
    
   case 'display'
    color.temp = [1 0 0]; % red for temporary orbits
    color.orb = [0 0 1];  % blue for orbits
    color.arrows = .4*[1 1 1]; % gray for arrows
    color.eqpt = [1 0 0];  % red for eq. pts.
    color.narrows = .4*[1 1 1];  % gray for nullcline arrows
    color.tx = [0 0 1]; % blue for xt plots & 3D plots
    color.ty = [1 0 0]; % red for yt plots
    color.ta = [1 0 0]; % red for axis plots
    color.sep = 2*[.5 0 .5];% [0,1,0];% green for separatrices
    color.xcline = [1 0 .75]; % magenta for xclines
    color.ycline = [1 .5 0]; %  orange for yclines
    color.level = 0.8*[.9,.5,.8];
    
   case 'bw'
    color.temp = [1 1 1]; % white for everything
    color.orb = [1 1 1];  
    color.eqpt = [1 1 1]; 
    color.arrows = [1 1 1];
    color.narrows = [1 1 1];  
    color.tx = [1 1 1]; 
    color.ty = [1 1 1]; 
    color.ta = [1 1 1]; 
    color.sep = [1 1 1];
    color.xcline = [1 1 1];
    color.ycline = [1 1 1];
    color.level = [1,1,1];
  end
  ud.color = color;
  ud.level = ' ';
  ppset = figure('name','pplane7 Setup','numb','off',...
		 'tag','pplane7','visible','off',...
		 'user',ud);
  
  pplane7('figdefault',ppset);
  frame(1) = uicontrol('style','frame','visible','off');
  eq(1)=uicontrol('style','text',...
		  'horizon','center',...
		  'string','The differential equations.','visible','off');
  ext = get(eq(1),'extent');
  rr=ext(4)/10;
  
  texth =ext(4)+4;    % 19;   		% Height of text boxes.
  varw = 40*rr;		% Length of variable boxes.
  equalw =13*rr;		% Length of equals.(30)
  eqlength = 230*rr;		% Length of right hand sides of equations.
  winstrlen = 120*rr;	% Length of string boxes in display frame.
  left = 0;		% Left margin of the frames.
  frsep = 1;    %3;		% Separation between frames.
  separation = texth;	% Separation between boxes.
  
  dfigwidth =2*left + varw+equalw+eqlength+10;	% Width of the figure.
  dfigurebot = 30;	% Bottom of the figure.
  buttw = dfigwidth/3;
  qwind = [0,frsep,buttw,separation];	% Quit button
  rwind = [buttw,frsep,buttw,separation];	% Revert "
  pwind = [2*buttw,frsep,buttw,separation];	% Proceed "
  
  disfrbot = 2*frsep + separation;	% Display frame.
  disfrw = winstrlen + varw +10;
  disfrht = 5*separation + 10;