Imperfect bifurcations

In this notebook, we will look into situations where the earlier symmetric cases, will include a parameter which will allow for imperfections, i.e. symmetry breaking.

In [10]:
import numpy as np;
import matplotlib.pyplot as plt;
plt.rcParams.update({"text.usetex":True});
%config InlineBackend.figure_format = "svg"
from ipywidgets import interactive
from scipy.optimize import fsolve

In the example below, we have two parameters $h$ and $r$, the former driving the breakage in symmetry.

In [2]:
def show_imperfect(h = 0.1, r = 0.5):

    def fx(x, h, r):
        return h + r*x - x**3

    x = np.linspace(-3, 3, 100);
    y = fx(x, h, r)
    plt.plot(x, y);
    plt.axhline(0, color='k');
    plt.ylim(-4, 4)
    plt.xlabel("x"); plt.ylabel("$\\dot{x}$")
w = interactive(show_imperfect, h = (-2, 2, 0.05), r = (-2, 2, 0.05))
w
In [3]:
def show_imp1(h, r):
    x = np.linspace(-3, 3, 100);
    y1 = x**3;
    y2 = h+r*x;
    plt.plot(x, y1, x, y2);
    plt.axhline(0, color='k');
    plt.ylim(-0.5, 0.5)
    
w = interactive(show_imp1, h = (-2, 2, 0.05), r = (-2, 2, 0.05))
w

Let us now draw the parameter space of (h, r).

In [11]:
from scipy.optimize import newton
In [12]:
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(-0.5, 0.5, 100);
h_a = np.linspace(-0.5, 0.5, 100);
[R, H] = np.meshgrid(r_a, h_a)
C = np.zeros(np.shape(R));

for i in np.arange(0, np.size(r_a)):
    for j in np.arange(0, np.size(h_a)):
        r = r_a[i];  h = h_a[j];
        xg = np.linspace(-3, 3, 5)
        
        sol = newton(f, fprime=dfdx ,x0=xg,  args=(r,h), full_output=True)
        roots = sol[0]
        whetherconv = sol[1]
        roots_conv = np.unique(np.round(roots[whetherconv], 3));
        number_roots = np.size(roots_conv);
        C[j,i] = number_roots;
    
plt.pcolor(R, H, C);    
ax = plt.gca(); ax.set_aspect(1);


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

F:\anaconda\lib\site-packages\scipy\optimize\zeros.py:459: RuntimeWarning: some failed to converge after 50 iterations
  warnings.warn(msg, RuntimeWarning)
In [13]:
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
In [14]:
from scipy.integrate import solve_ivp
In [15]:
def plot_imperfect(r=2, h = 0.24, x0 = 1.7, t_end =30):
    
    def xd(t, x, r, h):
        return h + r*x - x**3;
    
    sol = solve_ivp(xd, [0, t_end], [x0], args=[r,h], rtol=1e-6);
    plt.plot(sol.t, sol.y.T);
    
    plt.xlabel("$t$"); plt.ylabel("x");
    plt.title("$x = %f$"%(sol.y.T[-1]))

    
w = interactive(plot_imperfect,  r = (-0.5, 2, 0.01), h = (-2, 2, 0.01), x0 = (-2,2,0.1), t_end=(10, 50, 10));
w
In [9]:
from mayavi import mlab
h = np.linspace(-1,1);
r = np.linspace(-1,1);
x = np.linspace(-3,3)
mlab.plot3d(h, 0*h, 0*h, color=(1,0,0));
mlab.plot3d(0*r, r, 0*r, color=(0,1,0));
mlab.plot3d(0*x,0*x, x, color=(0,0,1));

H,R,X = np.meshgrid(h, r, x)
F = H + R*X - X**3;
G = X;

mlab.contour3d(H,R,X,F,contours=[-1,1])

src = mlab.pipeline.scalar_field(F)
src2 = mlab.pipeline.scalar_field(G)
mlab.pipeline.iso_surface(src, opacity=0.6,contours=[0],color=(1,0.4,0))
mlab.pipeline.iso_surface(src2,opacity=0.3,contours=[0])

mlab.axes(xlabel="R", ylabel="H", zlabel="X");
mlab.outline(True)
mlab.show()

No valid current object, please select an active object.
In [ ]: