Phase portraits in nonlinear systems

• In the previous lectures we have looked at behaviour of linear system and their phase plots and dependance of stability on the eigen values

• We now look at nonlinear systems

• Let us analyze the following system: $$\dot{x} = x + e^{-y}$$ $$\dot{y} = -y$$

• We reuse the code from prev lecture and plot the flow and trajectories

x =linspace(-3, 3, 20); y = linspace(-3, 3, 20);
[X, Y] = meshgrid(x, y);
U = X + exp(-Y);
V =-Y;
quiver(X, Y, U, V);
daspect([1 1 1]);
ylim([-3, 3]);
hold on;

x1 = linspace(-3, 3, 100); y1 = linspace(-3, 3, 100);
[X1, Y1] = meshgrid(x1, y1);
U1 = X1 + exp(-Y1);
V1 =-Y1;

xs = linspace(min(x1), max(x1), 20);
ys = linspace(min(y1), max(y1), 20);
[Xs, Ys] = meshgrid(xs,ys);

streamline(X, Y, U, V, [Xs], [Ys]);
contour(X1, Y1, U1, [0 0],color="g","linewidth",2);
contour(X1, Y1, V1, [0 0],color="m","linewidth",2);
hold off;

• Let us now consider the following system of non-linear equations: $$ \dot{x} = -x + x^3$$ $$ \dot{y} = -2y$$
• The fixed points are (0,0), (-1,0),(1,0)
• We can analyze the behaviour near the fixed points by linearizing the equations ,forming the Jacobian and finding the eigen values and eigen vectors

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=mysys(x, t) # returns the RHS
  y= ([-x(1) + x(1).^3,  -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 = -1.5+ 3*get (h.slider1, "value");
      y0 = -1.5 + 3*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+X.^3;
      V =-2*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 = -X1+X1.^3;
      V1 =-2*Y1;
      
      xs = linspace(min(x1), max(x1), 20);
      ys = linspace(min(y1), max(y1), 20);
      [Xs, Ys] = meshgrid(xs,ys);
      
      streamline(X1, Y1, U1, V1, [Xs], [Ys]);
      
      tspan=linspace(0,0.7,100)';    
      ics = [x0, y0];
      f = @(x,t) mysys(x, t);
      sol = lsode(f, ics, tspan);
      
      tout = linspace(0, max(tspan), 100);
      xout = sol(:,1);
      yout = sol(:,2);
      
      plot(xout,yout,"r","linewidth",2);
      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);

• Let us consider the problem: $$ \dot{x} = -y + ax(x^2+y^2)$$ $$ \dot{y} = x + ay(x^2 + y^2)$$
• The fixed point of the system is (0,0). When a=0 the trajectory is a closed path.
• When a>0, the trajectory spirals out and when a<0, the trajectory spirals in towards the fixed point.
• But if we linearize this system, then the eigen vectors indicate that the fixed point should have closed orbits around it irrespective of the value of a.
• So this case shows us that we cannot always linearize the equation to get the full picture.

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=mysys(x, t, a) # returns the RHS
  y= ([-x(2) + a*x(1).*(x(1).^2 + x(2).^2),  x(1) + a*x(2)*(x(1).^2 + x(2).^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;
    case {h.slider3}
      recalc = true;
    end
    
    if (recalc || init)  
      x0 = -1.5+ 3*get (h.slider1, "value");
      y0 = -1.5 + 3*get (h.slider2, "value");
      a = -0.4 + 0.8*get (h.slider3, "value");
      set (h.label1, "string", sprintf ("x0: %.1f", x0));
      set (h.label2, "string", sprintf ("y0: %.1f", y0));
      set (h.label3, "string", sprintf ("a: %1f", a));
      
      x = linspace(-3, 3, 20); y = linspace(-3, 3, 20);
      [X, Y] = meshgrid(x, y);
      U = -Y+a*X.*(X.^2+Y.^2);
      V = X + a*Y.*(X.^2 + Y.^2);
      quiver(X, Y, U./(U.^2+V.^2).^0.5, V./(U.^2+V.^2).^0.5);
      daspect([1 1 1]);
      hold on;
      
      tspan=linspace(0,4,100)';    
      ics = [x0, y0];
      f = @(x,t) mysys(x, t,a);
      sol = lsode(f, ics, tspan);
      
      tout = linspace(0, max(tspan), 100);
      xout = sol(:,1);
      yout = sol(:,2);
      
      plot(xout,yout,"r","linewidth",2);
      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]);
  
    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", 0.1,
  "position", [0.6 0.20 0.35 0.06]);
  
  set (gcf, "color", get(0, "defaultuicontrolbackgroundcolor"))
  guidata (gcf, h)
  update_plot (gcf, true);

• We now look at a conservative nonlinear system.
• We take a case of a oscillator with a non linear forcing function $$ m\ddot{x} = F $$ $$ F = x-x^3$$
• Let $y = \dot{x}$, now we have the following system of equations: $$ \dot{x} = y$$ $$ \dot{y} = x-x^3$$
• For this system we have 3 fixed points (0,0),(-1,0) and (1,0)
• We can analyze the behaviour near fixed points by making use of Jacobians

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=mysys(x, t) # returns the RHS
  y= ([x(2), x(1)-x(1).^3]);
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 = -1.5+ 3*get (h.slider1, "value");
      y0 = -1.5 + 3*get (h.slider2, "value");
      set (h.label1, "string", sprintf ("x0: %.1f", x0));
      set (h.label2, "string", sprintf ("y0: %.1f", y0));
      
      x = linspace(-3, 3, 15); y = linspace(-3, 3, 15);
      [X, Y] = meshgrid(x, y);
      U = Y;
      V = X - X.^3;
      quiver(X, Y, U./(U.^2+V.^2).^0.5, V./(U.^2+V.^2).^0.5);
      daspect([1 1 1]);
      hold on;
      
      tspan=linspace(0,20,100)';    
      ics = [x0, y0];
      f = @(x,t) mysys(x, t);
      sol = lsode(f, ics, tspan);
      
      tout = linspace(0, max(tspan), 100);
      xout = sol(:,1);
      yout = sol(:,2);
      
      plot(xout,yout,"r","linewidth",2);
      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);

• The energy of the system remains constant

function y=mysys(x, t) # returns the RHS
  y= ([x(2), x(1)-x(1).^3]);
end
tspan=linspace(0,20,100)';    
x0 = -1.1; y0 = 2.5;
ics = [x0, y0];
f = @(x,t) mysys(x, t);
sol = lsode(f, ics, tspan);

tout = linspace(0, max(tspan), 100);
xout = sol(:,1);
yout = sol(:,2);
E = 1/2*yout.^2 - 1/2*xout.^2 + 1/4*xout.^4;

#plt.plot(xout, yout)
plot(tout, E); 
ylim([1, 3]);

• We now vary the initial y coordinate and while keeping x coodinate equal to 1
• We plot in 3D by taking the z axis as the energy of the contour that the particle follows
  • We color code the orbits according to energy level

1;
function y=mysys(x, t) # returns the RHS
  y= ([x(2), x(1)-x(1).^3]);
end

tspan=linspace(0,20,100)';   
ys=linspace(0.1, 2.5, 30);
for y0=ys
  x0 = -1.1;
  ics = [x0, y0];
  f = @(x,t) mysys(x, t);
  sol = lsode(f, ics, tspan);
  
  tout = linspace(0, max(tspan), 100);
  xout = sol(:,1);
  yout = sol(:,2);
  E = 1/2*yout.^2 - 1/2*xout.^2 + 1/4*xout.^4;
  disp(mean(E))
  plot3(xout, yout, E);
  hold on;
end   
hold off;

-0.23398
-0.22227
-0.20373
-0.17833
-0.14608
-0.10698
-0.061038
-0.0082436
 0.051400
 0.11789
 0.19123
 0.27142
 0.35846
 0.45235
 0.55309
 0.66067
 0.77511
 0.89639
 1.0245
 1.1595
 1.3013
 1.4500
 1.6056
 1.7679
 1.9372
 2.1132
 2.2962
 2.4859
 2.6826
 2.8860