Chapter 8, Section 8.8.4, Figure 8.14a–c

Define the function h to be the Hamiltonian for the mixed resonance system (see Eq.(8.186)).

h[x_, y_, δ_, c_]:=
   (1/2)*δ*(x^2+y^2)+(1/4)*(x^2+y^2)^2+
   2*x*Sqrt[x^2+y^2+c]

Define the functions f and g as given in Eq.(8.185).

f[δ_,X_]:=δ+X^2;
g[i_, X_, c_]:=
(-1)^(i+1)*2*(Sqrt[X^2+c]/X + X/Sqrt[X^2+c]);

Calculate the derivatives of f and g.

df =  D[f[δ, X], X]

dg  =  D[g[0, X, c], X]

Plot both derivatives (for c = 1) as a function of X = R.

Plot[{df, (dg  /.  c -> 1)}, {X, 0, 10}]

Find the point of intersection.  At this point the tangents to f and g have equal gradient.

Xb  =  
      X  /.  FindRoot[
          df == (dg  /.  c -> 1),  {X, 1}]

Find the corresponding values of g and δb.

gval  =  g[0, Xb, 1]

δbval  =  gval  -  Xb^2

Plot f for \(TraditionalForm\`δ\&_\) = –15, –10, δbval, 0, 5 and 10, and g for c = 1, i = 0, 1.  This shows the tangential touching for \(TraditionalForm\`δ\&_\) =δbval, the bifurcation value (see Fig.8.14a).

\(p = 

      Plot[{f[\(-15\), X], f[\(-10\), X], f[δbval, \ X], \ f[0, X], 
          f[5, X], f[10, X], g[0, X, 1], g[1, X, 1]}, {X, 0, 4}, 
        PlotRange -> {{0, 4}, {\(-15\), 15}}, AxesOrigin -> {0, \(-15\)}, 
        PlotPoints -> 100, AspectRatio -> 1, \ Frame -> True]\)

Define a function, plotit, that uses Eq.(8.187) to plot the relationship between \(TraditionalForm\`δ\&_\) andthe value of x associated with the equilibrium points for the specified value c.  The function also plots a horizontal line for the \(TraditionalForm\`δ\&_\) associated with the bifurcation value.

\(plotit[c_, \ dbv_] := 

      Module[{}, 

        p1 = 
Plot[\((\(-x^3\) - 2*Sqrt[x^2 + c] - 
                  2*x^2/Sqrt[x^2 + c])\)/
              x, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {x, \(-4\), 
              4}, \ 
AspectRatio -> 1, \ 
            Frame -> True, 
PlotRange -> {{\(-4\), 
                  4}, {\(-15\), 15}}, DisplayFunction -> Identity]; 
        p2\  = \ Graphics[Line[{{\(-4\), dbv}, {4, dbv}}]]; 

        Show[p1, p2, DisplayFunction -> $DisplayFunction]]\)

Plot the function for c = 1 (see Fig.8.14b).

plotit[1, δbval]

Calculate the value of the Hamiltonian associated with the bifurcation point.

hval0 = h[Xb, 0, δbval,  1]

Plot the curves of constant Hamiltonian for h = –20, –10, –5, hval0 and 5 (see Fig.8.14c).

\(cp = 
      ContourPlot[
        h[x, y, δbval, \ 1], {x, \(-6\), 4}, {y, \(-5\), 
          5}, 
\ \ \ ContourShading -> False, 
        AxesOrigin -> {0, 0}, 
\ \ \ Contours -> {\(-20\), \(-10\), \(-5\), 
            hval0, 5}, AspectRatio -> 1, 
PlotPoints -> 150]\)

Plot h as a function of x for y = 0 and c = 1.

Plot[h[x, 0, δbval, 1], {x, -6, 4}]