Phase Transitions in Coupled Nonlinear Oscillators
Tanya Leise
Amherst College

Seek approximate solution x(t)=rcos(wt+f) 
for single oscillator x''-x'(a-bx'^2-gx^2)+w^2x=0
(exact solutions have period 2p/w +O(a^2))

ClearAll[x, alpha, beta, gamma, v, a, phi, phidot, r, rdot, 
      omega]
<< Graphics`

DE = 
      a -  v (alpha - beta  v^2 - gamma x^2 ) + 
        omega^2 x

x = r Cos[omega t + phi]
v = -r omega Sin[omega t + phi]
a = -rdot omega Sin[omega t + phi] - 
        r omega^2 Cos[omega t + phi] - 
        r omega phidot Cos[omega t + phi]

dots = 
      Simplify[Solve[{DE == 0, 
            rdot Cos[omega t + phi] - 
                phidot r Sin[omega t + phi] == 0}, {rdot, 
            phidot}]]

(* Average rdot over one period,  
      since assuming little short - 
        term change *) 
 (* 
      Including [[1]] removes curly braces from expression *) 
rdotaverage = Simplify[
          Integrate[
            rdot /. dots, {t, 0, 2 Pi/omega}]][[1]]

phidotaverage = Simplify[
          Integrate[
            phidot /. dots, {t, 0, 2 Pi/omega}]][[1]]

Solve[
      rdotaverage == 0, r]

Simplify[D[rdotaverage, r] /. r -> 0]

Simplify[
      D[rdotaverage, r] /. 
        r -> ((2 (Sqrt[alpha]))/(Sqrt[gamma + 3 beta (omega^2)]))]

(* Since d (rdot)/dr is positive when r = 0,  
      this solution is unstable. *) 
 (* On the other hand,  
      d (rdot)/
            dr is negative when r = ((2 (Sqrt[alpha]))/(Sqrt[gamma + 3 beta (omega^2)])),  so this solution is stable. *)

alpha = 0.05;
beta = 0.007;
gamma = 12;
omega = 1;

xnSoln = 
      NDSolve[{xn''[
                t] -  xn'[
                  t] (alpha - beta xn'[t]^2 - gamma xn[t]^2 ) + 
              omega^2 xn[t] == 0, 
          xn[0] ==  .05, xn'[0] == 0}, 
        xn[t], {t, 0, 100 Pi}]

Plot[{2  Sqrt[alpha]/Sqrt[gamma + 3 beta omega^2] Cos[omega t], 
        xn[t] /. xnSoln}, {t, 0 Pi, 50 Pi}, 
      PlotStyle -> {Hue[0], Hue[ .7]}, 
      PlotLegend -> {"Periodic Approximation", "Numerical solution\
"}, LegendSize -> {1.9,  .2}, LegendPosition -> {-1, -1}, 
      AxesLabel -> {"t", "x"}]

Seek approximate solutions for coupled oscillators
Assume all parameters positive with alpha>2a and gamma>8b.

ClearAll[xL, xR, alpha, beta, gamma, vL, vR, aL, aR, phiL, phiR, 
      phiLdot, phiRdot, rL, rR, rLdot, rRdot, omega, a, b]

deL = 
      aL -  vL (alpha - beta vL^2 - gamma xL^2) + 
        omega^2 xL + (vR - vL) (-a + b (xR - xL)^2); 
    deR = aR -  vR (alpha - beta  vR^2 - gamma xR^2) + 
        omega^2 xR - (vR - vL) (-a + 
              b (xR - xL)^2);

xL = rL Cos[omega t + phiL];
xR = rR Cos[omega t + phiR];
vL = -rL omega Sin[
            omega t + phiL];
vR = -rR omega Sin[
            omega t + phiR];
aL = -omega Cos[omega t + phiL] rL (omega + phiLdot) - 
        omega Sin[omega t + phiL] rLdot
aR = -omega Cos[omega t + phiR] rR (omega + phiRdot) - 
        omega Sin[omega t + phiR] rRdot

DEL = Simplify[deL]
DER = Simplify[deR]

consistencyR = -rR Sin[omega t + phiR] phiRdot + 
        Cos[omega t + phiR] rRdot
consistencyL = -rL Sin[omega t + phiL] phiLdot + 
        Cos[omega t + phiL] rLdot

dots = 
      Simplify[Solve[{DEL == 0, DER == 0, 
            consistencyL == 0, consistencyR == 0}, {rLdot, rRdot, 
            phiLdot, phiRdot}]]

rLdotaverage = Simplify[
          Integrate[
            rLdot /. dots, {t, 0, 
              2 Pi/omega}]][[1]]

    rRdotaverage = Simplify[
          Integrate[
            rRdot /. dots, {t, 0, 
              2 Pi/omega}]][[1]]

    phiLdotaverage = Simplify[
          Integrate[
            phiLdot /. dots, {t, 0, 
              2 Pi/omega}]][[1]]

    phiRdotaverage = Simplify[
          Integrate[
            phiRdot /. dots, {t, 0, 2 Pi/omega}]][[1]]

rLdotave = 
      rLdotaverage /. {phiL -> phi + phiR}

    rRdotave = rRdotaverage /. {phiL -> phi + phiR}

phidotave = 
      Simplify[phiLdotaverage - phiRdotaverage] /. {phiL -> 
            phi + phiR}

Solve[(rLdotave /. {phi -> 0, 
              rL -> r, rR -> r}) == 0, r]

Solve[(rLdotave /. {phi -> 
                Pi, rL -> r, rR -> r}) == 0, 
      r]

Solve[(rLdotave /. {phi -> 
                ArcCos[((-2 a + b (r^2))/(b (r^2)))], rL -> r, 
              rR -> r}) == 0, r]

J = {{Simplify[
              D[rLdotave, rL]], Simplify[D[rLdotave, rR]], 
            Simplify[D[rLdotave, phi]]}, {Simplify[D[rRdotave, rL]], 
            Simplify[D[rRdotave, rR]], 
            Simplify[D[rRdotave, phi]]}, {Simplify[D[phidotave, rL]], 
            Simplify[D[phidotave, rR]], 
            Simplify[D[phidotave, phi]]}};

(* 
      In - phase oscillatory solution is stable since all eigenvalues are negative *) 
Jinphase = 
        Simplify[
          J /. {phi -> 0, 
              rL -> ((2 (Sqrt[alpha]))/(Sqrt[gamma + 3 beta (omega^2)])), 
              rR -> ((2 (Sqrt[alpha]))/(Sqrt[gamma + 3 beta (omega^2)]))}];

    MatrixForm[Jinphase]

    Eigenvalues[Jinphase]

(* 
      Out - of - 
        phase oscillatory solution may or may not be stable,  
      depending on omega *) 
Joutphase = 
        Simplify[
          J /. {phi -> Pi, 
              rL -> ((2 (Sqrt[alpha - 2  a]))/(Sqrt[-8 b + gamma + 3 beta (omega^2)])), 
              rR -> ((2 (Sqrt[alpha - 2  a]))/(Sqrt[-8 b + gamma + 3 beta (omega^2)]))}];

    MatrixForm[Joutphase]

    Simplify[
      Eigenvalues[{{Joutphase[[1, 1]], 
            Joutphase[[1, 2]]}, {Joutphase[[2, 1]], 
            Joutphase[[2, 2]]}}]]

(* 
      Eigenvalue d (phidot)/
            d (phi) = ((2 (-4 alpha b + a (gamma + 3 beta (omega^2))) π)/(omega (-8 b + gamma + 3 beta (omega^2)))) can change signs,  
      while other two eigenvalues are always negative. *) 
omegacritical = 
      Solve[((2 (-4 alpha b + a (gamma + 3 beta (omega^2))) π)/(omega (-8 b + gamma + 3 beta (omega^2)))) == 0, 
        omega]

(* Check stability of third possibility;  
      note double angle formula used to simplify trig expressions *) 
 (* 
      When this solution exists (omega < 
            omegacritical) it is unstable since there is a positive eigenvalue *) 
Jother = 
        Simplify[J /. {Cos[2 phi] -> 
                    2  Cos[phi]^2 - 1} /. {phi -> 
                  ArcCos[((-4 a + b (rL^2) + b (rR^2))/(2 b rL rR))]} /. {rL -> ((2 (Sqrt[alpha]))/(Sqrt[gamma + 3 beta (omega^2)])), 
              rR -> ((2 (Sqrt[alpha]))/(Sqrt[gamma + 3 beta (omega^2)]))}];

    MatrixForm[Jother]

    Eigenvalues[Jother]