dfield7.m

matlab of linear system

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

% dfield7	is an interactive tool for studying single first order
%	differential equations.  When dfield7 is executed, a dfield7 Setup
%	window is opened.  The user may enter the differential
%	equation and specify a display window using the interactive
%	controls in the Setup window.
%
%	When the Proceed button is pressed on the Setup window, the DF
%	Display window is opened.  At first this window displays a
%	direction line field for the differential equation.  When the
%	mouse button is depressed in the dfield7 Display window, the
%	solution to the differential equation with that initial
%	condition is calculated and plotted.
%
%	Other options are available in the Options menu.  These are
%	fairly self explanatory.  The Settings option allows the user
%	to change several parameters.  Included here are the
%	possibilities of using a vector field instead of the default
%	line field, and of changing the number of field points computed
%	and displayed.

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


  startstr = 'dfield7';

  if nargin < 1
    action ='initialize';
  end


  % disp(action)


  if strcmp(action,'initialize')
    
    % First we make sure that there is no other copy of dfield7
    % running, since this causes problems.	
    
    dfh = findobj('tag','dfield7');
    if ~isempty(dfh);
      qstring = {'There are some open dfield7 figures (perhaps invisible).  ';...
		 'What do you want to do?'};
      tstring = 'Only one copy of dfield7 can be open at one time.';
      b1str = 'Restart dfield7.';
      b2str = 'Just close those figures.';
      b3str = 'Do nothing.';
      answer = questdlg(qstring,tstring,b1str,b2str,b3str,b1str);
      if strcmp(answer,b1str)
	delete(dfh);
	eval(startstr);return
      elseif strcmp(answer,b2str)
	delete(dfh);return
      else
	return
      end 
    end  % if ~isempty(dfh);
    
    % 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 dfield7.
    
    oldfiles = dir('dftp*.m');
    kk = zeros(0,1);
    for k = 1:length(oldfiles)
      fn = oldfiles(k).name;
      fid = fopen(fn,'r');
      if fid ~= -1
	ll = fgetl(fid);
	ll = fgetl(fid);
	ll = fgetl(fid);
	fclose(fid);
	if strcmp(ll,'%% Created by dfield7')
	  delete(fn)
	else
	  kk = [kk;k];
	end
      end
    end
    oldfiles = dir([tempdir,'dftp*.m']);
    for k = 1:length(oldfiles)
      fn = [tempdir,oldfiles(k).name];
      fid = fopen(fn,'r');
      if fid ~= -1
	ll = fgetl(fid);
	ll = fgetl(fid);
	ll = fgetl(fid);
	fclose(fid);
	if isempty(findstr('%% Created by DF',ll))
	  kk = [kk;k];
	else
	  delete(fn)
	end
      end
    end
    if ~isempty(kk)
      fprintf('The following files\n');
      for j = 1:length(kk)
	fn = [tempdir,oldfiles(kk(j)).name];
	fn = strrep(fn,'\','\\');
	fprintf(['     ',fn,'\n']);
      end
      fprintf('may have been created by an older version of DFIELD.\n');
      fprintf('If so they should be deleted.\n\n');
    end
    
    style = 'white';
    dfdir = pwd;
    ssize = get(0,'defaultaxesfontsize');
    npts = 20;
    solver = 'Dormand-Prince';
    solvhandle = @dfdp45;
    if exist('dfstart','file')
      H = dfstart;
      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,'dfdir')
	  dfdir = H.dfdir;  
	end
	if isfield(H,'solver')
	  solver = H.solver;  
	  switch solver
	   case 'Dormand-Prince'
	    solvhandle = @dfdp45;
	    
	   case 'Euler'	
	    solvhandle = @dfeul;
	    
	   case 'Runge-Kutta 2'	
	    solvhandle = @dfrk2;
	    
	   case 'Runge-Kutta 4'	
	    solvhandle = @dfrk4;
	    
	   otherwise
	    error('Undefined solver.');
	  end
	end
      end
    end
    if get(0,'screendepth') == 1
      style = 'bw';
    end
    ud.ssize = ssize;
    ud.dfdir = dfdir;
    comp = computer;
    if strcmp(comp,'PCWIN')
      ud.fontsize = 0.8*ud.ssize;
    else
      ud.fontsize = ud.ssize;
    end
    system.name = 'default equation';
    system.xname = 'x';
    system.tname = 't';
    system.der = ' x^2 - t';
    system.wind = [-2 10 -4 4];
    system.pname = {};
    system.pval = {};
    
    system(2).name = 'logistic equation';
    system(2).xname = 'P';
    system(2).tname = 't';
    system(2).der = ' r*P*(1 - P/K)';
    system(2).wind = [0 20 0 15];
    system(2).pname = {'r','K','',''};
    system(2).pval = {'0.75','10','',''};
    
    system(3).name = 'RC circuit';
    system(3).xname = 'V_c';
    system(3).tname = 't';
    system(3).der = '(A*cos(\omega*t) - V_c)/(R*C)';
    system(3).wind = [0 10 -5 5];
    system(3).pname = {'A'  '\omega'  'R'  'C'};
    system(3).pval = {'5'  '3'  '0.5'  '2'};

    system(4).name = 'RL circuit';
    system(4).xname = 'I';
    system(4).tname = 't';
    system(4).der = '(A*cos(\omega*t) - R*I)/L';
    system(4).wind = [0 10 -5 5];
    system(4).pname = {'A'  '\omega'  'R'  'L'};
    system(4).pval = {'5'  '3'  '3'  '2'};

    ud.c = system(1);	% Changed values.
    ud.o = system(1);	% Original values.
    ud.fieldtype = 'lines';
    ud.npts = npts;  
    ud.style = style;
    ud.settings.magn = 1.25;
    ud.settings.refine = 4;
    ud.settings.tol = 1e-6;
    ud.settings.solver = solver;
    ud.settings.solvhandle = solvhandle;
    ud.settings.stepsize = 0.1;
    ud.settings.hmax = 0;
    ud.settings.speed = 100;
    ud.system = system;
    switch style
     case 'black'
      color.temp = [1 0 0]; % red for temporary orbits
      color.orb = [1 1 0];  % yellow for orbits
      color.arrows = .5*[1 1 1];  % gray for arrows
      color.pline = [1 1 1];    
      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.arrows = .7*[1 1 1]; % gray for arrows
      color.pline = [0 0 0];    
      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.pline = [0 0 0];    
      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.pline = [0 0 0];   
      color.level = [255 150 10]/255;%  0.8*[.9,.5,.8];
      
     case 'bw'
      color.temp = [1 1 1]; % white for everything
      color.orb = [1 1 1];  
      color.arrows = [1 1 1];
      color.pline = [1 1 1];    
      color.level = [1,1,1];
    end
    ud.color = color;
    
    dfset = figure('name','dfield7 Setup','numb','off',...
		   'tag','dfield7','visible','off',...
		   'user',ud);
    
    dfield7('figdefault',dfset);
    frame(1) = uicontrol('style','frame','visible','off');
    eq(1)=uicontrol('style','text',...
		    'horizon','center',...
		    'string','The differential equation.','visible','off');
    ext = get(eq(1),'extent');
    rr=ext(4)/10;
    
    texth =ext(4)+4;      % 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 = 1;		% Left margin of the frames.
    frsep = 1;    	% Separation between frames.
    separation = texth;	% Separation between boxes.
    
    disfrw = varw+equalw+eqlength+10;
    
    dfigwidth =2*left + disfrw;	% Width of the figure.
    dfigurebot = 60;	% 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.
    disfrht = 3*separation + 10;
    disfrwind = [left, disfrbot, disfrw, disfrht];
    
    pfrbot = disfrbot + disfrht + frsep;
    pfrw = disfrw;
    pfrht = 2*separation + 10;
    pfrwind = [left, pfrbot, pfrw, pfrht];
    
    defrbot = pfrbot + pfrht + frsep;	% Equation frame.
    defrw = pfrw;
    defrht = 3*separation + 15;
    defrwind = [left, defrbot, defrw, defrht];
    
    dfigureheight = defrbot + defrht +frsep;	% Height of the figure.
    
    set(dfset,'pos',[30 dfigurebot dfigwidth dfigureheight]);
    
    set(frame(1),'pos',defrwind);
    
    
    xname=[
	'ud = get(gcf,''user'');'...
	'Xname=get(ud.h.xname,''string'');'...
	'set(ud.h.twind(3),''string'','...
	'[''The minimum value of '',Xname,'' = '']);'...
	'set(ud.h.twind(4),''string'','...
	'[''The maximum value of '',Xname,'' = '']);'...
	'ud.c.xname = Xname;'...
	'ud.flag = 0;'...
	'set(gcf,''user'',ud);'];
    
    tname=[
	'ud = get(gcf,''user'');'...
	'Tname=get(ud.h.tname,''string'');'...
	'set(ud.h.twind(1),''string'','...
	'[''The minimum value of '',Tname,'' = '']);'...
	'set(ud.h.twind(2),''string'','...
	'[''The maximum value of '',Tname,'' = '']);'...
	'ud.c.tname = Tname;'...
	'ud.flag = 0;'...
	'set(gcf,''user'',ud);'];
    
    der =[
	'ud = get(gcf,''user'');'...
	'ud.c.der = get(ud.h.der,''string'');'...
	'ud.flag = 0;'...
	'set(gcf,''user'',ud);'];
    
    equationbot = defrbot + 5;
    eqlabelbot = equationbot + 2*separation;
    derbot = equationbot + separation + 5;		% Bottom of equation.
    tbot = equationbot;				% Bottom of ind. var.
    
    lablen =200*rr;
    
    eqlableft = (dfigwidth-lablen)/2;
    eqleft = left + 5;
    fudge = 0.15*separation;
    set(eq(1),'pos',[eqlableft eqlabelbot lablen texth]);
    tcolor = get(gcf,'defaultuicontrolbackgroundcolor');
    ecolor = 'w';
    
    ud.h.xname=uicontrol('pos',[eqleft, derbot, varw, texth],...
			 'style','edit',...
			 'horizon','right',...
			 'string',ud.o.xname,...
			 'call',xname,...
			 'visible','off',...
			 'backgroundcolor',ecolor);
    
    eq(2) = uicontrol('style','text',...
		      'pos',[eqleft+varw derbot-fudge equalw texth],...
		      'horizon','center','string',''' = ','visible','off',...
		      'backgroundcolor',tcolor);
    
    ud.h.der=uicontrol('pos',[eqleft+varw+equalw derbot eqlength texth],...
		       'string',ud.o.der,...
		       'horizon','left','style','edit',...
		       'call',der,'visible','off','backgroundcolor',ecolor);
    
    eq(3) = uicontrol('style','text',...
		      'horizon','right',...
		      'string','The independent variable is ',...
		      'visible','off','backgroundcolor',tcolor);
    
    ext = get(eq(3),'extent');
    tnw = ext(3)+10;
    tvarl = eqleft + tnw;
    set(eq(3),'pos',[eqleft,tbot-fudge,tnw,texth]);
    
    ud.h.tname = uicontrol('pos',[tvarl tbot varw texth],...
			   'style','edit',...
			   'horizon','left',...
			   'string',ud.o.tname,...
			   'call',tname,...
			   'visible','off',...
			   'backgroundcolor',ecolor);
  
    pframe = uicontrol('style','frame','pos',pfrwind,'visible','on');
    
    pncall = [
	'[h,fig] = gcbo;'...
	'ud =  get(fig,''user'');'...
	'num = get(h,''user'');'...
	'ud.c.pname{num} = get(ud.h.pname(num),''string'');'...
	'ud.flag = 0;'...
	'set(gcf,''user'',ud);'];
    
    pvcall = [
	'[h,fig] = gcbo;'...
	'ud =  get(fig,''user'');'...
	'num = get(h,''user'');'...
	'ud.c.pval{num} = get(ud.h.pval(num),''string'');'...
	'ud.flag = 0;'...
	'set(gcf,''user'',ud);'];
    
    pnamew = 30*rr;
    peqw = 8*rr;
    pbot1 = pfrbot + 5;
    pbot2 = pbot1 +separation;
    pbot(1) = pbot2;
    pbot(2) = pbot1;
    
    
    paratit=uicontrol('style','text',...
		      'horizon','center',...
		      'string',{'Parameters';'&';'expressions:'},...
		      'visible','off','backgroundcolor',tcolor);
    
    ext = get(paratit,'extent');
    paratitw = ext(3);
    pos = [eqleft pbot(2) paratitw 2*texth];
    set(paratit,'pos',pos);
    
    psep = 20;
    pvalw = (dfigwidth - 2*eqleft - paratitw)/2 - psep - pnamew - peqw;
    pval = ud.c.pval;
    pname = ud.c.pname;
    for jj = 1:2
      for kk = 1:2
	pleft = eqleft + paratitw +psep +(kk-1)*(pnamew+peqw+pvalw+ psep);
	peqleft = pleft + pnamew;
	pvleft = peqleft + peqw;
	K = kk +2*(jj-1);
	if K > length(pname)
	  pname{K} = '';
	end
	name = pname{K};
	if K <= length(pval)
	  value = pval{K};
	else
	  value = '';
	end
	ud.h.pname(K) = uicontrol('style','edit',...
				  'pos',[pleft pbot(jj) pnamew texth],...
				  'horizon','right','string',name,...
				  'user',K,...
				  'call',pncall,...
				  'visible','off',...
				  'backgroundcolor','w');
	equal(K) = uicontrol('style','text',...
			     'pos',[peqleft pbot(jj)-fudge peqw texth],...
			     'horizon','center',...