Processing math: 100%

Linear 2D flows contd

In general, for any matrix relating $\mathbf{x}$ to $\dot{\mathbf{x}}$, say $\mathbf{A}$, we can find out the eigenvalues and eigenvectors. Based on that we can analyze various trajectories in the phase space and understand the role of the eigenproperties of the matrix.

clear;
clc; close all;

h.ax = axes ("position", [0.05 0.4 0.5 0.5]);  #reduce plot windows size to accomodate sliders

function update_plot (obj, init = false)
  ## gcbo holds the handle of the control
  clc;
  h = guidata (obj);
  replot = false;
  recalc = false;
  switch (gcbo)                       # If we make any change then we replot
    case {h.print_pushbutton}
      fn =  uiputfile ("*.png");
      print (fn);
    case {h.slider1}
      recalc = true;
    case {h.slider2}
      recalc = true;
     case {h.slider3}
      recalc = true;
     case {h.slider4}
      recalc = true;
    end
    
    if (recalc || init)  
      a = -3 + 6*get (h.slider1, "value");
      b = -3 + 6*get (h.slider2, "value");
      c =  -3 + 7*get (h.slider3, "value");
      d = -3 + 6*get(h.slider4, "value");
      set (h.label1, "string", sprintf ("a: %1f", a));
      set (h.label2, "string", sprintf ("b: %1f", b));
      set (h.label3, "string", sprintf ("c: %1f", c)); 
      set (h.label4, "string", sprintf ("d: %1f", d));
      
    #a = 1; b = 1; c = 4; d = -2;
    A = [[a b]; [c d]];
    [v, lam] = eig(A);
    disp(lam);
    disp(v(:,1)');
    disp(v(:,2)');
    

    x = linspace(-3, 3, 20); y = linspace(-3, 3, 20);
    [X, Y] = meshgrid(x, y);
    U = a*X + b*Y;
    V = c*X + d*Y;
    quiver(X, Y, U, V);
    daspect([1 1 1]);
    hold on;
    x1 = linspace(-3, 3, 100); y1 = linspace(-3, 3, 100);
    [X1, Y1] = meshgrid(x1, y1);
    U1 = a*X1 + b*Y1;
    V1 = c*X1 + d*Y1;
    
    seedpoints = [x1 x1;zeros(size(y1))+3 zeros(size(y1))-3];
    streamline(X1, Y1, U1, V1,seedpoints(1,:),seedpoints(2,:));
    ylim([-3, 3]);
    plot(x1, x1.*(v(2,1)./v(1,1)));
    plot(x1, x1.*(v(2,2)./v(1,2)));
    hold off;
      
      guidata (obj, h);
    else
      set (h.plot, "ydata", y);
    end
  end
  
  # specify the GUI plot properties like sliders and label formatting
  
  ## print figure
  h.print_pushbutton = uicontrol ("style", "pushbutton",
  "units", "normalized",
  "string", "print plot\n(pushbutton)",
  "callback", @update_plot,
  "position", [0.6 0.45 0.35 0.09]);
  ## guess
  h.label1 = uicontrol ("style", "text",
  "units", "normalized",
  "string", "Guess:",
  "horizontalalignment", "left",
  "position", [0.05 0.25 0.35 0.08]);
  
  h.slider1 = uicontrol ("style", "slider",
  "units", "normalized",
  "string", "slider",
  "callback", @update_plot,
  "value", 1,
  "position", [0.05 0.20 0.35 0.06]);
  
  h.label2 = uicontrol ("style", "text",
  "units", "normalized",
  "string", "Iteration:",
  "horizontalalignment", "left",
  "position", [0.05 0.15 0.35 0.08]);
  
  h.slider2 = uicontrol ("style", "slider",
  "units", "normalized",
  "string", "slider",
  "callback", @update_plot,
  "value", 1,
  "position", [0.05 0.10 0.35 0.06]);
  
    h.label3 = uicontrol ("style", "text",
  "units", "normalized",
  "string", "Guess:",
  "horizontalalignment", "left",
  "position", [0.6 0.25 0.35 0.08]);
  
  h.slider3 = uicontrol ("style", "slider",
  "units", "normalized",
  "string", "slider",
  "callback", @update_plot,
  "value", 1,
  "position", [0.6 0.20 0.35 0.06]);
  
  h.label4 = uicontrol ("style", "text",
  "units", "normalized",
  "string", "Iteration:",
  "horizontalalignment", "left",
  "position", [0.6 0.15 0.35 0.08]);
  
  h.slider4 = uicontrol ("style", "slider",
  "units", "normalized",
  "string", "slider",
  "callback", @update_plot,
  "value", 1,
  "position", [0.6 0.10 0.35 0.06]);
  
  set (gcf, "color", get(0, "defaultuicontrolbackgroundcolor"))
  guidata (gcf, h)
  update_plot (gcf, true);

• We see that for complex eigen values, with positive real part, the values spiral out and the magnitude grows.
• For real part negative, the values spiral in.
• When the real part is zero, we have a closed orbit
• When both eigen values become equal, all the vectors move out from the origin. They either grow exponentially or decay depending upon the sign of eigen value.

• Let us now plot the trajectory for a=1, b=1, c=-2, d=44
• We make use of lsode routine to solve the initial value problem
  • We provide the initial condition and the time span

clear; clc; close all;

h.ax = axes ("position", [0.05 0.4 0.5 0.5]);  #reduce plot windows size to accomodate sliders
function y = f(x, t)
  y =  [x(1) + x(2), 4*x(1) - 2*x(2)];
end  


function update_plot (obj, init = false)
  ## gcbo holds the handle of the control
  h = guidata (obj);
  replot = false;
  recalc = false;
  switch (gcbo)                       # If we make any change then we replot
    case {h.print_pushbutton}
      fn =  uiputfile ("*.png");
      print (fn);
    case {h.slider1}
      recalc = true;
    case {h.slider2}
      recalc = true;
    end
    
    if (recalc || init)  
      x0 = -2.5 + 5*get (h.slider1, "value");
      y0 = -2.5 + 5*get (h.slider2, "value");
      set (h.label1, "string", sprintf ("x0: %1f", x0));
      set (h.label2, "string", sprintf ("y0: %1f", y0));
      
      x = linspace(-3,3,20);
      y = linspace(-3,3,20);
      [X, Y] = meshgrid(x,y);
      U = X + Y;
      V = 4*X - 2*Y;
      quiver(X, Y, U, V);
      daspect([1 1 1]);
      hold on;
            
      tspan=linspace(0,1,100)';
      xd = @(x,t) f(x, t);
      
      sol = lsode(xd, [x0, y0], tspan);
      tout = linspace(0, max(tspan), 100);
      xout = sol(:,1);
      yout = sol(:,2);
      plot(xout, yout);
      xout
      yout
      plot(x, x, '--k');
      plot(x, -4*x, '-k');
      plot(x0, y0, 'ok');
      xlim([-3,3]);               
      ylim([-3, 3]);
      guidata (obj, h);
      hold off;
    else
      set (h.plot, "ydata", y);
    end
  end
  
  # specify the GUI plot properties like sliders and label formatting
  
  ## print figure
  h.print_pushbutton = uicontrol ("style", "pushbutton",
  "units", "normalized",
  "string", "print plot\n(pushbutton)",
  "callback", @update_plot,
  "position", [0.6 0.45 0.35 0.09]);
  ## guess
  h.label1 = uicontrol ("style", "text",
  "units", "normalized",
  "string", "Guess:",
  "horizontalalignment", "left",
  "position", [0.05 0.25 0.35 0.08]);
  
  h.slider1 = uicontrol ("style", "slider",
  "units", "normalized",
  "string", "slider",
  "callback", @update_plot,
  "value", 0.1,
  "position", [0.05 0.20 0.35 0.06]);
  
  h.label2 = uicontrol ("style", "text",
  "units", "normalized",
  "string", "Iteration:",
  "horizontalalignment", "left",
  "position", [0.05 0.15 0.35 0.08]);
  
  h.slider2 = uicontrol ("style", "slider",
  "units", "normalized",
  "string", "slider",
  "callback", @update_plot,
  "value", 0.5,
  "position", [0.05 0.10 0.35 0.06]);
  
  set (gcf, "color", get(0, "defaultuicontrolbackgroundcolor"))
  guidata (gcf, h)
  update_plot (gcf, true);