Una interface gráfica para estudiar las trayectorias en plano de fase de dos estados.

function dfase(op)

% DFASE(OP) Es una interface gráfica para estudiar las trayectorias en plano de fase
% de dos estados.
%
% Tiene la posibilidad de graficar:
%   => Las trayectoria de los estados ingresando la condición inicial con la
%   posición del mouse cuando se hace un click.
%   => Las isoclinas.
%   => Los vectores que indican la pendiente.
%
% Haciendo un click con el mouse en el plano de fase, se cálculan y grafican las
% trayectorias que describen los estados.
% El método de resolución está basado en el método de Runge-Kutta, consultar
% help de ODE45 para mayores destalles. Las ecuaciones pueden ser lineales o no.
% Como tiene una tolerancia fijada de antemano puede que algunas ecuaciones no lineales
% muy complicadas fallen en la convergencia. En la mayoría de los casos no presentan
% problemas.
%
% Para utilizarla simplemente ingrese "dfase3" desde el Command Window
%
%


% si no se ingresa opción entonces se define la GUI
if nargin<1
   op=1;
end,


if op==1
   % dibujo toda la GUI
 
   % para evitar problemas de localización de objetos
   h=findobj('Tag','FigFase');
if ~isempty(h)
  errordlg('Ya existe una ventana Diagrama de Fase',...
     'Error en Diagrama de Fase',...
     'modal');
  return,
end,

 
   Hf = figure('Units','normalized',...
    'Position',[0.1 0.1 0.8 0.8],...
  'NumberTitle','off',...
  'Name','Diagrama de Fase',...
      'Menubar','none',...
      'Tag','FigFase');
 
   c=get(gcf,'Color');
 
% plano de fase
Ha1=axes('Units','normalized',...
  'Position',[0.04 0.25 0.6 0.69],...
  'FontSize',8,...
  'Xlim',[-10 10],...
      'Ylim',[-10 10],...
      'box','on',...
  'Tag','fase',... % el tag es "fase"
  'NextPlot','add',...
  'ButtonDownFcn','dfase(3)');

title('Diagrama de fase');
 
   % gráfico de la velocidad (x2)
Ha2=axes('Units','normalized',...
'Position',[ 0.68 0.25 0.3 0.3 ],...
      'FontSize',8,...
      'box',' on',...
  'NextPlot','add',...
    'Tag','velocidad'); % el tag es "velocidad"

title('Velocidad (x2)');

% gráfico de la posición (x1)
Ha3=axes('Units','normalized',...
  'Position',[ 0.68 0.64 0.3 0.3 ],...
      'box',' on',...
  'FontSize',8,...
  'NextPlot','add',...
  'Tag','posicion');

title('Posición (x1)');
 
   % segunda ecuación
uicontrol('Style','text',...
      'Units','normalized',...
      'BackgroundColor',c,...
      'ForegroundColor',[1 1 1],...
  'Position',[ 0.05 0.15 0.2 0.03 ],...
  'HorizontalAlignment','left',...
  'String','Primer ecuación:')

uicontrol('Style','edit',...
  'Units','normalized',...
  'BackgroundColor',[1 1 1],...
  'Position',[ 0.25 0.14 0.2 0.05 ],...
  'HorizontalAlignment','left',...
    'String','x2',...
  'Tag','PrimerEc')

   % segunda ecuación
uicontrol('Style','text',...
  'Units','normalized',...
      'BackgroundColor',c,...
      'ForegroundColor',[1 1 1],...
  'Position',[ 0.05 0.08 0.2 0.03 ],...
  'HorizontalAlignment','left',...
  'String','Segunda ecuación:')
 
uicontrol('Style','edit',...
  'Units','normalized',...
  'BackgroundColor',[1 1 1],...
  'Position',[ 0.25 0.07 0.2 0.05 ],...
    'HorizontalAlignment','left',...
  'String','-x1-x2',...
  'Tag','SegundaEc');

% boton vectores
uicontrol('Style','push',...
  'Units','normalized',...
  'Position',[ 0.5 0.14 0.12 0.05 ],...
  'String','Vectores',...
  'Tag','Push1',...
  'Callback','dfase(2)');

% boton isoclinas
uicontrol('Style','push',...
  'Units','normalized',...
  'Position',[ 0.5 0.07 0.12 0.05 ],...
  'String','Isoclinas',...
  'Tag','Push2',...
  'Callback','dfase(5)');
 
% boton borrar
uicontrol('Style','push',...
  'Units','normalized',...
  'Position',[ 0.65 0.07 0.12 0.05 ],...
  'String','Borrar',...
  'Tag','Push3',...
  'Callback','dfase(4)');
 
   % boton imprimir
uicontrol('Style','push',...
  'Units','normalized',...
  'Position',[ 0.65 0.14 0.12 0.05 ],...
  'String','Imprimir',...
  'Tag','Push4',...
  'Callback','print');

   % boton cerrar
uicontrol('Style','push',...
  'Units','normalized',...
  'Position',[ 0.8 0.07 0.12 0.05 ],...
  'String','Cerrar',...
  'Tag','Push5',...
  'Callback','close(gcf)');

% grafico dos líneas que pasan por el origen
set(Hf,'CurrentAxes',Ha1);
plot([-10 0;10 0],[0 -10;0 10],...
  'Color',[1 1 1],...
  'LineWidth',0.1);

elseif op==2
   % Para graficar los vectores
 
   % obtengo las ecuaciones
   ec1=get(findobj('Tag','PrimerEc'),'String');
   ec2=get(findobj('Tag','SegundaEc'),'String');
 
   % planteo la ecuación de la pendiente
   ec=['m=(' ec2 ')./(' ec1 ');'];

   % hago al axes del plano de fase el actual
   Ha=findobj('Tag','fase');
   set(gcf,'CurrentAxes',Ha);
 
% grafico los vectores
for j=1:20
      for s=1:20
         % calculo la pendiente y para un dado punto del plano
         a=1.5;b=0.75;
         x1=-10+j*a;
     x2=-10+s*a;
         eval(ec);
         % calculo el punto final de modo que tenga una longitud fija
         xd=sqrt((b^2)/(1+m^2));
      x21=m*xd+x2;
      x11=x1+xd;
  plot([x1 x11],[x2 x21],'y')
end,
end,
 
 
elseif op==3
   % Para calcular y graficar las trayectorias de los estados
   % utilizo ODE45 (Runge-Kuta)
 
   % primero necesito armar un archivo .m con las ecuaciones diferenciales
   f=fopen('ectmp.m','wt+');
 
   fprintf(f,'%s','function yd=ecuacion(t,x)');
   fprintf(f,'\n\n\n\n');
 
   % acomodo las ecuaciones ingresadas
   ec1=get(findobj('Tag','PrimerEc'),'String');
   ec2=get(findobj('Tag','SegundaEc'),'String');
 
   % necesito exprezarlas como un vector para la ODE45
   ec11=strrep(ec1,'x1','x(1)');
   ec22=strrep(ec2,'x1','x(1)');
   ec11=strrep(ec11,'x2','x(2)');
   ec22=strrep(ec22,'x2','x(2)');

   ec=['yd=[' ec11 ';' ec22 '];'];
 
   fprintf(f,'%s\n',ec);
 
   fclose(f);
 
   % capturo las condiciones iniciales desde la posición actual del mouse
   yo=get(findobj('Tag','fase'),'CurrentPoint');
 
   % resuelvo las ecuaciones diferenciales
   [t,y]=ode45('ectmp',0,10,yo(2,1:2)');
 
   % elimino la función que tiene las ecuaciones de la memoria
   % para que los cálculos sean actuales
   clear ectmp
 
   % grafico
   set(gcf,'CurrentAxes',findobj('Tag','velocidad'));
   hl = findobj('Parent',gca,'Type','line');
   delete(hl);
   plot(t,y(:,1),'w');
 
   set(gcf,'CurrentAxes',findobj('Tag','posicion'));
   hl = findobj('Parent',gca,'Type','line');
   delete(hl);
   plot(t,y(:,2),'w');
 
   set(gcf,'CurrentAxes',findobj('Tag','fase'));
   plot(y(:,1),y(:,2),'w');
 
 
elseif op==4
   % Para borrar los vectores, isoclinas y trayectorias de estados
 
   % hago al axes del plano de fase el actual
   Ha=findobj('Tag','fase');
   set(gcf,'CurrentAxes',Ha);
 
   % busco y borro los objetos tipo líneas
   Hdel=findobj('Parent',Ha,'Type','line');
   delete(Hdel);
 
   % redibujo las líneas que pasan por el origen
   plot([-10 0;10 0],[0 -10;0 10],...
    'Color',[1 1 1],...
  'LineWidth',0.1);


elseif op==5
   % Para calcular y graficar las isoclinas
 
% obtengo las ecuaciones que se ingresaron en forma de cadena  
   ec1=get(findobj('Tag','PrimerEc'),'String');
   ec2=get(findobj('Tag','SegundaEc'),'String');
 
   % trabajo las ecuaciones para trabajar punto a punto
   ec=['m=(' ec2 ')/(' ec1 ')'];
   ec1=strrep(ec,'*','.*');
   ec1=strrep(ec1,'/','./');
   ec1=strrep(ec1,'^','.^');
 
   % calculo 20 valores de pendientes de manera que queden distribuidos
   % en el palno de fase
   x1=[-10*ones(1,5) -10 -6 -2 2 6 10*ones(1,5) 10 6 2 -2 -6];
   x2=[-10 -6 -2 2 6 10*ones(1,5) 10 6 2 -2 -6 -10*ones(1,5)];
   eval([ec1 ';']);
   M=m;
 
   % resuelvo la ecuación para poner en forma explícita x2
   % SE NECESITA SYMBOLIC TOOLBOX
   ecs=solve(ec,'x2');
ecs1=ecs(1,:);
   ecs1=strrep(ecs1,'*','.*');
   ecs1=strrep(ecs1,'/','./');
   ecs1=strrep(ecs1,'^','.^');
 
   % calculo numericamente las isoclinas
x1=-10:0.2:10;
for j=1:20;
      m=M(j);
eval(['x2=' ecs1 ';']);
plot(x1,x2,...
    'Color',[0.8 0.8 0.8])
end,
 
 
end