Bifurcations and 3D flows

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
  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;
    end
    
    if (recalc || init)  
      mu = -1 + 2*get (h.slider1, "value");
      set (h.label1, "string", sprintf ("mu: %1f", mu));
      x = linspace(-3, 3, 100); y = linspace(-3, 3, 100);
      [X, Y] = meshgrid(x, y);
      U = mu-X.^2;
      V = -Y;
      x1 = linspace(min(x), max(x), 20);
      y1 = linspace(min(y), max(y), 20);
      [X1, Y1] = meshgrid(x1,y1);
      streamline(X, Y, U, V, [X1], [Y1]);
      daspect([1 1 1]);
      hold on;
      line(xlim, [0 0],"color", "m");
      line([0 0], ylim,"color", "g");    
      plot(mu.^0.5, 0, 'ok',"markersize",10);
      plot(-mu.^0.5, 0, 'sr',"markersize",10);
      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]);
  
  set (gcf, "color", get(0, "defaultuicontrolbackgroundcolor"))
  guidata (gcf, h)
  update_plot (gcf, true);

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
  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)  
      a = -2 + 4*get (h.slider1, "value");
      b = -2 + 4*get (h.slider2, "value");
      set (h.label1, "string", sprintf ("a: %1f", a));
      set (h.label2, "string", sprintf ("b: %1f", b));
      x = linspace(-3, 3, 100); y = linspace(-3, 3, 100);
      [X, Y] = meshgrid(x, y);
      U = -a*X+Y;
      V = -X.^2/(1+X.^2) -b*Y;
      x1 = linspace(min(x), max(x), 20);
      y1 = linspace(min(y), max(y), 20);
      [X1, Y1] = meshgrid(x1,y1);
      streamline(X, Y, U, V, [X1], [Y1]);
      daspect([1 1 1]);
      hold on;
      line(xlim, [0 0],"color", "k");
      line([0 0], ylim,"color", "k");    
      contour(X, Y, U, [0 0],color="g","linewidth",2);
      contour(X, Y, V, [0 0],color="m","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.8,
  "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.8,
  "position", [0.05 0.10 0.35 0.06]);
  
  set (gcf, "color", get(0, "defaultuicontrolbackgroundcolor"))
  guidata (gcf, h)
  update_plot (gcf, true);

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
  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;
    end
    
    if (recalc || init)  
      mu = -1 + 2*get (h.slider1, "value");
      set (h.label1, "string", sprintf ("mu: %1f", mu));
      x = linspace(-3, 3, 30); y = linspace(-3, 3, 30);
      [X, Y] = meshgrid(x, y);
      U = mu*X-Y+X.*Y.^2;
      V = X+mu*Y+Y.^3;
      x1 = linspace(min(x), max(x), 20);
      y1 = linspace(min(y), max(y), 20);
      [X1, Y1] = meshgrid(x1,y1);
      streamline(X, Y, U, V, [X1], [Y1]);
      daspect([1 1 1]);
      hold on;
      line([0 0], ylim,"color","k"); line(xlim, [0 0],"color","k");
      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]);
  
  set (gcf, "color", get(0, "defaultuicontrolbackgroundcolor"))
  guidata (gcf, h)
  update_plot (gcf, true);

clc;
clear; close all;

function y=mysys(x, t, gamma, d, omega) # returns the RHS
    y=([x(2), x(1)-x(1).^3-gamma*x(2)+d*cos(omega*t)]);
end

gamma = 0.1; d = 0.35; omega = 1.4; tp = 2*pi/omega;
Nmax = 30000;
tmax = Nmax*tp;
x0 = 0; 
y0 = 0;
ics = [x0, y0];
tspan= linspace((Nmax-50)*tp, tmax, 2000);
f = @(x,t) mysys(x, t, gamma, d, omega);
sol = lsode(f, ics, tspan);

plot(sol(:,1), sol(:,2));

clc;
clear; close all;

function y=mysys(x, t, gamma, d, omega) # returns the RHS
    y=([x(2), x(1)-x(1).^3-gamma*x(2)+d*cos(omega*t)]);
end

gamma = 0.1; d = 0.35; omega = 1.4; tp = 2*pi/omega;
Nmax = 30000;
tmax = Nmax*tp;
x0 = 0; 
y0 = 0;
ics = [x0, y0];
tspan= tp*(Nmax-20000:1:Nmax);
f = @(x,t) mysys(x, t, gamma, d, omega);
sol = lsode(f, ics, tspan);

plot(sol(:,1), sol(:,2),'.k', "markersize",1);

** Warning: Cannot interrupt Octave on Windows

clear; clc;
function y=mysys(x, t, s, r, b) # returns the RHS
    y=([s*(x(2)-x(1)), r*x(1)-x(2)-x(1)*x(3), x(1)*x(2)-b*x(3)]);
end
tmax = 200;
tspan = linspace(0, tmax, 20000);

s = 10; b = 8/3; r = 28;
x0 = 0; 
y0 = 1;
z0 = 0;
ics = [x0, y0, z0];
f = @(x,t) mysys(x, t, s, r, b);
sol = lsode(f, ics, tspan);

xout = sol(:,1);
yout = sol(:,2);
zout = sol(:,3);
plot(xout, zout);

r = (xout.^2+yout.^2+zout.^2).^0.5;
patch(xout, yout, zout, r,'FaceColor','none','EdgeColor','interp');

clear; clc;

pkg load signal

function y=mysys(x, t, s, r, b) # returns the RHS
    y=([s*(x(2)-x(1)), r*x(1)-x(2)-x(1)*x(3), x(1)*x(2)-b*x(3)]);
end
tmax = 200;
tspan = linspace(0, tmax, 20000);

s = 10; b = 8/3; r = 28;
x0 = 0; 
y0 = 1;
z0 = 0;
ics = [x0, y0, z0];
f = @(x,t) mysys(x, t, s, r, b);
sol = lsode(f, ics, tspan);

xout = sol(:,1);
yout = sol(:,2);
zout = sol(:,3);

[peaks idx]= findpeaks(zout);

subplot(2,1,1);
plot(tspan, zout);
daspect([2 1]);
hold on;
plot(tspan(idx), peaks, 'bx');
z = peaks;
zk = z(1:end-2); zkp1 = z(2:end-1);
hold off;
subplot(2,1,2);
plot(zk, zkp1, '.k', "markersize",5);
hold on;
daspect([1 1 1]);
xd = linspace(30, 45);
plot(xd, xd, '-m');
hold off;

clear; clc;
function y=mysys(t, x, s, r, b) # returns the RHS
    y=([s*(x(2)-x(1)), r*x(1)-x(2)-x(1)*x(3), x(1)*x(2)-b*x(3)]);
end

function [position,isterminal,direction] = ycross(t, x, s, r, b)
    position =  x(2);
    isterminal = 0;  % Halt integration 
    direction = 0;   % The zero can be approached from either direction
end

tmax = 200;
tspan = linspace(0, tmax, 20000);

s = 10; b = 8/3; r = 28;
x0 = 0; 
y0 = 1;
z0 = 0;
ics = [x0, y0, z0];
f = @(t,x) mysys(t, x, s, r, b);
funcEvent = @(t,x) ycross(t, x, s, r, b);
opt = odeset("Events", funcEvent,"RelTol",1e-6); 
sol = ode45(f, tspan, ics, opt);

xout = sol.y(1,:);
yout = sol.y(2,:);
zout = sol.y(3,:);
plot(xout, zout);

te = sol.xe;
xe = sol.ye(1,:);
ye = sol.ye(2,:);
ze = sol.ye(3,:);
r =  1.0*(yout>0.0);
patch(xout, yout, zout, r,'FaceColor','none','EdgeColor','interp');
hold on;
scatter3(xe, ye, ze,2,[1 1 1]);
hold off;