Imperfect bifurcations

• We had previously seen the behavior of $\dot{x} = rx-x^3$
• The parameter k can break the symmetry of the above equation
  • $\dot{x} = k + rx - x^3$

• We define a function tak takes as input the parameters k and r
• We then plot the function $\dot{x}$ w.r.t x

set(0, "defaultlinelinewidth", 2);
set (0, "defaulttextfontname", "TimesNewRoman")
set (0, "defaulttextfontsize", 20)
set (0, "DefaultAxesFontName", "TimesNewRoman")
set(0, 'DefaultAxesFontSize', 20)

close all
clear h

graphics_toolkit qt

h.ax = axes ("position", [0.05 0.4 0.5 0.5]); #reduce plot window 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)
    case {h.print_pushbutton}
      fn =  uiputfile ("*.png");
      print (fn);
    case {h.slider1}
      recalc = true;
    case {h.slider2}
      recalc = true;
  end
  
  if (recalc || init)
    k = -2+ 4*get (h.slider1, "value"); # change the range of frequency slider
    r = -2 + 4*get (h.slider2, "value"); # change the range of points slider
    x = linspace(-3, 3, 100);
    set (h.label1, "string", sprintf ("k: %.1f", k));
    set (h.label2, "string", sprintf ("r: %.1f", r));      
    y = k + r*x - x.^3;
      
    h.plot = plot(x, y, x, 0.*x, '-k');
    h.ylim= ylim([-4 4]);
    guidata (obj, h);
    else
      set (h.plot, "ydata", y);
   end
end


## 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]);
h.label1 = uicontrol ("style", "text",
                           "units", "normalized",
                           "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.5,
                            "position", [0.05 0.20 0.35 0.06]);
                            
h.label2 = uicontrol ("style", "text",
                           "units", "normalized",
                           "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);

• It is not exactly clear how r and h determine the number of roots
• Let us split the equation into two functions $h + rx$ and $x^3$
• We now look at the intersection of these two curves to find the fixed points

set(0, "defaultlinelinewidth", 2);
set (0, "defaulttextfontname", "TimesNewRoman")
set (0, "defaulttextfontsize", 20)
set (0, "DefaultAxesFontName", "TimesNewRoman")
set(0, 'DefaultAxesFontSize', 20)

close all
clear h

graphics_toolkit qt

h.ax = axes ("position", [0.05 0.4 0.5 0.5]); #reduce plot window 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)
    case {h.print_pushbutton}
      fn =  uiputfile ("*.png");
      print (fn);
    case {h.slider1}
      recalc = true;
    case {h.slider2}
      recalc = true;
  end
  
  if (recalc || init)
    k = -2+ 4*get (h.slider1, "value"); # change the range of frequency slider
    r = -2 + 4*get (h.slider2, "value"); # change the range of points slider
    x = linspace(-3, 3, 100);
    set (h.label1, "string", sprintf ("k: %.1f", k));
    set (h.label2, "string", sprintf ("r: %.1f", r));
    y1 = x.^3;
    y2 = k+r*x;
    h.plot = plot(x, y1, x, y2, x, 0.*x, '-k');      
    h.ylim= ylim([-0.5 0.5]);
    guidata (obj, h);
    else
      set (h.plot, "ydata", y);
   end
end


## 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]);
h.label1 = uicontrol ("style", "text",
                           "units", "normalized",
                           "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.5,
                            "position", [0.05 0.20 0.35 0.06]);
                            
h.label2 = uicontrol ("style", "text",
                           "units", "normalized",
                           "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 now draw the parameter space of (h, r) and identify regions according to number of roots
• We color code the diagram according to the number of roots
• We also use newton solver instead of fsolve
• We store the unique roots returned by newton solver and then plot the stability diagram

function y=f(x, r, h)
    y=h + r*x- x.^3;
end

r_a = linspace(-0.5, 0.5, 100);
h_a = linspace(-0.5, 0.5, 100);
[R, H] = meshgrid(r_a, h_a);
C = zeros(size(R));

for i=1:size(r_a)(2)
    for j=1:size(h_a)(2)
        r = r_a(i);  h = h_a(j);
        xg = linspace(-3, 3, 5); 
    disp(i)
    disp(j)
    f = f(x,r,h);
        [roots, res] = fsolve(f, xg);
        roots_conv = unique(round(roots(whetherconv<1e-4), 3));
        number_roots = size(roots_conv);
        C(j,i) = number_roots;
    end
end
  
pcolor(R, H, C);  shading interp;  
ax = gca(); ax.set_aspect(1);


xlabel("r");
ylabel("h");
colorbar();
title("Stability diagram");

 1
 1
error: 'x' undefined near line 7 column 11
error: scalar cannot be indexed with .

• Cusp shapes are highlights of imperfect bifurcations

• Now we look at the stability of the roots
• We plot the bifurcation diagram wrt to r and we vary h using the widget

def show_hvar(h):
    def f(x, r, h):
        return h + r*x- x**3;

    def dfdx(x, r, h):
        return r-3*x**2;

    r_a = np.linspace(-2, 2,20);


    for r in r_a:
        sol = newton(f, fprime=dfdx ,x0=xg,  args=(r,h), full_output=True)
        for i in np.arange(0, np.size(sol[0])):
            if sol[1][i] == True: # Solution has convered
                root = sol[0][i];
                slope_at_root = dfdx(root, r, h)
                if slope_at_root > 0:
                    plt.plot(r, root, 'bx')
                else:
                    plt.plot(r, root, 'ro')

    ax = plt.gca(); ax.set_aspect(1);
    plt.xlim(np.min(r_a), np.max(r_a))
    plt.axhline(0); plt.axvline(0)
    plt.xlabel("r");
    plt.ylabel("x*");
    plt.title("Imperfect bifurcation for: h = %3.2f"%h);
    r1 = np.linspace(-2, 2);
    x1 = np.linspace(-2, 2);
    [R1,X1] = np.meshgrid(r1, x1);
    F  = h+ R1*X1 - X1**3;

    plt.contour(R1, X1, F, levels=[0.0])
    


w = interactive(show_hvar, h = (-2, 2, 0.01))
w

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, r, k)
  y =  k + r*x - x.^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;
     case {h.slider3}
      recalc = true;
     case {h.slider4}
      recalc = true;
    end
    
    if (recalc || init)  
      x0 = -2 + 4*get (h.slider1, "value");
      r = -0.5 + 2.5*get (h.slider2, "value");
      t_end =  10 + 40*get (h.slider3, "value");
      k = -2 + 4*get(slider4, "value");
      set (h.label1, "string", sprintf ("x0: %.1f", x0));
      set (h.label2, "string", sprintf ("r: %.1f", r));
      set (h.label3, "string", sprintf ("t_end: %.1f", t_end)); 
      set (h.label4, "string", sprintf ("k: %.1f", k));
      t=linspace(0,t_end,100)';
      xd = @(x,t) f(x, t, r,k);
      sol = lsode(xd, x0, t);
      plot(t, sol, "-k");
      xlabel("t"); ylabel("x");               
      ylim([-1.2, 1.2]);
      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", 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]);
  
  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", 0.5,
  "position", [0.6 0.10 0.35 0.06]);
  
  set (gcf, "color", get(0, "defaultuicontrolbackgroundcolor"))
  guidata (gcf, h)
  update_plot (gcf, true);