The Damped Driven Pendulum

Jesse Welch

Plan of Concentration 
Marlboro College 2007

Initialization Functions

These functions are each called later by other functions and thus are required no matter which section is evaluated beyond here.

Setting the parameter values and derivative functions for the damped driven pendulum.

Code which defines the function

SetPendParam[driverforce_, driverfreq_, dampening_, driverphase_] := 
        Module[{}, 
dθ[t_, x_, y_, z_] = 
            y; 
dω[t_, x_, y_, z_] = 
            driverforce * Cos[z + driverphase] - Sin[x] - 
              dampening * y; 
dϕ[t_, x_, y_, z_] = 
            driverfreq; 
frequency  = 
            driverfreq; 
F[
              t_, {x_, y_, z_}] = {dθ[t, x, y, z], 
              dω[t, x, y, z], dϕ[t, x, y, z]};];

Explaination and definitions

This function sets the parameter values for the differential equations to be solved.  There are four parameter values which define the pendulum system and three variable time-dependant  functions which are our solutions.  There is no direct output when this function is run.
	Local Input Values:
					driveforce - value which sets the amplitude of the driving force function
					drivefreq - value which sets the frequency of the driving force function
					dampening -  value which sets the linear dampening constant 
					driverphase -   value which sets the initial phase of the driving force function
	Global Parameter Defined:
					dθ - First order derivative equation for θ[t],  a linear equation in terms of ω
					dω - First order derivative equation for ω[t], a combination of nonlinear terms of θ and ϕ, as well as a negative linear term for ω
					dϕ - First order derivative equation for ϕ[t],  a constant
					frequency - constant set equal to the drive frequency input value, called later to set length of period
					F - function of t and the set {θ, ω, ϕ}, which creates a list of the three derivative values numerically evaluated as a function of t, θ, ω, and ϕ

Examples

As an example we set the parameters {driveforce = 0.9, drivefreq = 2/3, dampening = 0.5, drivephase = 0} :

SetPendParam[0.9, 2/3, 0.5, 0]

We have now set the function F[t, {θ, ω, ϕ}], which we can call to see how it was defined :

F[t, {θ, ω, ϕ}]

Integrating the equations of motion using a fourth order Runge-Kutta method giving only one point per cycle.

Code which defines the function

RungeOncePerCycle[initialt_, initialθ_, initialω_, 
          stepspercycle_, 
          drivecycles_] := 
Module[{j, 
            stepsper = stepspercycle, 
            cycles = drivecycles}, 
period = 
            2.0*N[π, 
                20]*(1/frequency);  (*                        defines the period of the drive cycles *) 
h = (period/stepsper);  (*                                                              defines the stepsize *) 
trajectory  =  
            Table[{}, {cycles}];  (*             Empty table which will be filled by the algorithm *) 
T = 
            Table[(initialt + (j - 1)* h), {j, 1, 
                stepsper}];  (*                  Table defining the time values*) 
R = {initialθ, 
              initialω, 
              initialt};  (*                                     List of initial conditions*) 
Do[{ (*                                               Loop which takes the solution forward one step*) 
For[j = 1, 
                j < stepsper + 1, j++, 
k1 = 
                  h*F[k*T[[j]], R]; 
k2 = 
                  h*F[k*(T[[j]] + h/2), 
                      R + k1/2]; 
k3 = 
                  h*F[k*(T[[j]] + h/2), 
                      R + k2/2]; 
k4 = 
                  h*F[k*(T[[j]] + h), 
                      R + k3]; 
Rnew = 
                  R + 1/6 (k1 + 2 k2 + 2 k3 + k4);  (* 
                  Average the four values for the slopes and define the next R*) 
R = 
                  Rnew;
]; 
trajectory[[k]] = 
                R;  (*             Add current value of the three variables to the solution table*) 
}, {k, 1, 
              cycles}]
];

Explaination and definitions

This function is a Runge-Kutta fourth order differential equation solving algorithm.  It uses the global function F defined from the SetPendParam function and uses it's values to solve for the variable functions given in the equations.  
	Local Input Values:
					initialt - value which sets the initial value of t to solve from
					initialθ - value which sets the initial value of θ to solve from
					initialω -  value which sets the initial value of ω to solve from
	Hidden Local Values:
	  				stepspercycle -   value which sets the number of steps to take per drivecycle ( where drivecycles are a measure of elapsed time)
					divecycles - value which sets the total number of drivecycles to integrate forward in time 
	Global Parameter Called:	
					frequency - is called on from SetPendParam
					F	     - is called on from SetPendParam
	Global Parameter Defined:
					period - defined to be the value of the length of a drivecycle, calls on frequency
					h - defined to be the length of a single step, defined by period divided by stepspercycle
					trajectory - table created initially as only 0's, filled in during the algorithm by the solution values to the three dependant variables {θ, ω, ϕ}
					T - table created to contain the value of time at each step during a drivecycle
					R - list of initial values of {θ, ω, ϕ} given by input values defined above, and then set equal to proceeding point, to use again as initial values
					k1 - defined to be the slope of the solutions at the initial point of the time step using the value of y at the initial point
					k2 - defined to be the slope of the solutions at the halfway point of the time step using the value of k1 as the initial point
					k3 - defined to be the slope of the solutions at the halfway point of the time step using the value of k2 as the initial point
					k4 - defined to be the slope of the solutions at the end point of the time step using k3 as the initial point
					Rnew - defines the value of the solution at the end point of a time step, taking a weighted average of the four k's.  
					
This function saves only one data point per drivecycle, at the end of the last timestep of each drivecycle, as some calculations done require saving only one piece of data per cycle.  There is no direct output, though a global list trajectory is defined, these sets are the solutions to the differential equation, given as a list of three piece vectors in the form {θ, ω, ϕ}.

Examples

As an example we call RungeOncePerCycle with initial values {θ, ω, ϕ} = {0, 0, 0}, 100 steps per drivecycle, over a range of 10 drivecycles:

RungeOncePerCycle[0, 0, 0, 100, 10]

Note: This only works here because we have already run the SetPendParam function with some parameter values to describe the differential equation system.

We then can see the results as a table of 10 lists, each list with 3 values, each of the three values represents the value of {θ, ω, ϕ} at the end of each drivecycle, and we integrated for 10 drivecycles, we end with 10 sets in our final table.

trajectory

Integrating the equations of motion using a fourth order Runge-Kutta method giving all points in the cycle.

Code which defines the function

RungeAll[initialt_, initialθ_, initialω_, 
          stepspercycle_, 
          drivecycles_] := 
Module[{j, 
            stepsper = stepspercycle, 
            cycles = drivecycles}, 
period = 
            2.0*N[π, 20]*
              frequency;   
h = 
            (period/stepsper);   
trajectory  =  
            Table[{}, {stepsper}];   
trajectory2  =  
            Table[{}, {cycles}];   
T = 
            Table[(initialt + (j - 1)* h), {j, 1, 
                stepsper}];   
R = {initialθ, 
              initialω, 
              initialt};   
Do[{  
For[j = 1, 
                j < stepsper + 1, j++, 
k1 = 
                  h*F[k*T[[j]], R]; 
k2 = 
                  h*F[k*(T[[j]] + h/2), 
                      R + k1/2]; 
k3 = 
                  h*F[k*(T[[j]] + h/2), 
                      R + k2/2]; 
k4 = 
                  h*F[k*(T[[j]] + h), 
                      R + k3]; 
Rnew = 
                  R + 1/6 (k1 + 2 k2 + 2 k3 + 
                          k4);   
R = 
                  Rnew; 
trajectory[[j]] = 
                  
                  R;   
]; 
trajectory2[[k]] = 
                trajectory;   
}, {k, 1, cycles}]; 
Do[
            trajectory2[[g + 1]] = 
              Join[trajectory2[[g]], trajectory2[[g + 1]]], {g, 1, 
              cycles - 1}]; 
trajectory2  =  Drop[
                trajectory2, cycles - 1][[1]];   
];

Explaination and definitions

This function is a Runge-Kutta fourth order differential equation solving algorithm.  It uses the global function F defined from the SetPendParam function and uses it's values to solve for the variable functions given in the equations.  
	Local Input Values:
					initialt - value which sets the initial value of t to solve from
					initialθ - value which sets the initial value of θ to solve from
					initialω -  value which sets the initial value of ω to solve from
	Hidden Local Values:
	  				stepspercycle -   value which sets the number of steps to take per drivecycle ( where drivecycles are a measure of elapsed time)
					divecycles - value which sets the total number of drivecycles to integrate forward in time 
	Global Parameter Called:	
					frequency - is called on from SetPendParam
					F	     - is called on from SetPendParam
	Global Parameter Defined:
					period - defined to be the value of the length of a drivecycle, calls on frequency
					h - defined to be the length of a single step, defined by period divided by stepspercycle
					trajectory - table created initially as only 0's, filled in during the algorithm `by the solution values to the three dependant variables {θ, ω, ϕ} at each 							point through each cycle, each cycle overwrites the previous
					trajectory2 - table created initially as only 0's, filled in during the algorithm by the table of one drivecycle solution with values of the three 								dependant variables {θ, ω, ϕ}, doing this for each cycle we end up with the entire set of cycles in one table 
					T - table created to contain the value of time at each step during a drivecycle
					R - list of initial values of {θ, ω, ϕ} given by input values defined above, and then set equal to proceeding point, to use again as initial values
					k1 - defined to be the slope of the solutions at the initial point of the time step using the value of y at the initial point
					k2 - defined to be the slope of the solutions at the halfway point of the time step using the value of k1 as the initial point
					k3 - defined to be the slope of the solutions at the halfway point of the time step using the value of k2 as the initial point
					k4 - defined to be the slope of the solutions at the end point of the time step using k3 as the initial point
					Rnew - defines the value of the solution at the end point of a time step, taking a weighted average of the four k's.  
					
This function saves all data points per drivecycle, as some calculations done require saving every piece of data per cycle.  There is no direct output, though a global list trajectory2 is defined, these sets are the solutions to the differential equation, given as a list of three piece vectors in the form {θ, ω, ϕ}.

Examples

As an example we call RungeAll with initial values {θ, ω, ϕ} = {0, 0, 0}, 100 steps per drivecycle, over a range of 10 drivecycles:

RungeAll[0, 0, 0, 100, 10]

We could then see the results as a table of 1000 lists, each list with 3 values, each of the three values represents the value of {θ, ω, ϕ} at each point of each drivecycle, and we integrated for 10 drivecycles, we end with 1000 sets in our final table, the whole thing wil  not be shown.

{trajectory2[[100]], trajectory2[[468]]}

Plotting Functions

These functions output graphical plots in one form or another...

Plotting Theta vs t

Code which defines the function

PlotTheta[driverforce_, dampening_, stepspercycle_, drivecycles_, 
        numofsolutions_] := 
Module[{stepsper = 
            stepspercycle, 
          cycles = drivecycles}, 
SetPendParam[driverforce, 
          2/3, dampening, 
          0]; 
trajec3 = {{0, 
              0}}; 
Do[{
a = 
              Random[Real, {0, 2  π}]; 
aa = 
              Random[Real, {0, 
                  2  π}]; 
trajectory3  =  
              Table[{}, {stepsper * cycles}]; 
RungeAll[
              0, a, aa, stepsper, cycles]; 
Do[
              trajectory3[[k]]  =  {trajectory2[[k]][[3]], 
                  Mod[trajectory2[[k]][[1]] - π, 
                    2  π]}
, {k, 1, 
                stepsper * cycles}]; 
trajec3 = 
              Join[trajec3, 
                trajectory3];}, {numofsolutions}]; 
ListPlot[trajec3];
]

Explaination and definitions

This function plots the solution to the set of differential equations defined by parameter values input into this function in terms of θ vs t.  Note: the two parameters, drivephase and drivefreq, are set equal to, 0 and 2/3, respectively in all plots to reduce the number of inputs required to plot said function, also initial conditions for θ and ω are randomly selected for the same reasons.  
	Local Input Values:
					driverforce - value which sets the amplitude of the driving force function
					dampening -  value which sets the linear dampening constant 
					numofsolutions - value which sets the number of solutions to calculate for this given parameter set
	Hidden Local Values:
					stepspercycle -   value which sets the number of steps to take per drivecycle ( where drivecycles are a measure of elapsed time)
					divecycles - value which sets the total number of drivecycles to integrate forward in time 
	Global Parameters Called: 
					there are no direct references to any global parameters previously defined 
	Global Functions Called:	
					SetPendParam - called to set the equations with parameter values input into this function
					RungeAll - called to calculate the values of {θ, ω, ϕ} at each point over the input integration time
					ListPlot - called to output the table of data computed by RungeAll
	Global Parameter Defined:
					a - random real value in the range {0, 2π}which sets the initial value of θ
					aa - random real value in the range {0, 2π}which sets the initial value of ω
					trajec3 - table which contains values for all combined solutions
					trajectory3 - table which contains values for one solution, then written over to contain the second solution, etc.

Examples

As an example, we can call this function with driverforce = 1.081, dampening = 0.5, stepspercycle = 150, cycles = 40, numofsolutions = 1, and get a listplot of one solution of θ vs t over 40 drivecycles, for a total of 6000 total points evaluated, and thus we aquire a plot the solution for the parameter set {1.081, 0.5, 2/3, 0}

PlotTheta[1.081, 0.5, 150, 30, 1]

Plotting Omega vs t

Code which defines the function

PlotOmega[driverforce_, dampening_, stepspercycle_, drivecycles_, 
        numofsolutions_] := 
Module[{stepsper = 
            stepspercycle, 
          cycles = drivecycles}, 
γ  =  
          numofsolutions; 
SetPendParam[driverforce, 2/3, 
          dampening, 
          0]; 
trajec3 = {{0, 
              0}}; 
Do[{
a = 
              Random[Real, {0, 2  π}]; 
aa = 
              Random[Real, {0, 
                  2  π}]; 
trajectory3  =  
              Table[{}, {stepsper * cycles}]; 
RungeAll[
              0, a, aa, stepsper, cycles]; 
Do[
              trajectory3[[k]]  =  {trajectory2[[k]][[3]], trajectory2[[k]][[2]]}, {k, 1, 
                stepsper * cycles}]; 
trajec3 = 
              Join[trajec3, 
                trajectory3];}, {γ}]; 
ListPlot[
          trajec3];
]

Explaination and definitions

This function plots the solution to the set of differential equations defined by parameter values input into this function in terms of ω vs t.  Note: the two parameters, drivephase and drivefreq, are set equal to, 0 and 2/3, respectively in all plots to reduce the number of inputs required to plot said function, also initial conditions for θ and ω are randomly selected for the same reasons.  
	Local Input Values:
					driverforce - value which sets the amplitude of the driving force function
					dampening -  value which sets the linear dampening constant 
					numofsolutions - value which sets the number of solutions to calculate for this given parameter set
	Hidden Local Values:
					stepspercycle -   value which sets the number of steps to take per drivecycle ( where drivecycles are a measure of elapsed time)
					divecycles - value which sets the total number of drivecycles to integrate forward in time 
	Global Parameters Called: 
					there are no direct references to any global parameters previously defined 
	Global Functions Called:	
					SetPendParam - called to set the equations with parameter values input into this function
					RungeAll - called to calculate the values of {θ, ω, ϕ} at each point over the input integration time
					ListPlot - called to output the table of data computed by RungeAll
	Global Parameter Defined:
					a - random real value in the range {0, 2π}which sets the initial value of θ
					aa - random real value in the range {0, 2π}which sets the initial value of ω
					trajec3 - table which contains values for all combined solutions
					trajectory3 - table which contains values for one solution, then written over to contain the second solution, etc.

Examples

As an example, we can call this function with driverforce = 1.081, dampening = 0.5, stepspercycle = 150, cycles = 40, numofsolutions = 1, and get a listplot of one solution of ω vs t over 40 drivecycles, for a total of 6000 total points evaluated, and thus we aquire a plot the solution for the parameter set {1.081, 0.5, 2/3, 0}

PlotOmega[1.081, 0.5, 150, 40, 1]

Plotting Theta vs Omega (Phase Plots)

Code which defines the function

PlotPhase[driverforce_, dampening_, stepspercycle_, drivecycles_, 
        numofsolutions_, 
        killtransients_] := 
Module[{stepsper = 
            stepspercycle, 
          cycles = drivecycles}, 
γ  =  
          numofsolutions; 
SetPendParam[driverforce, 2/3, 
          dampening, 
          0]; 
trajec3 = {{0, 
              0}}; 
Do[{
a = 
              Random[Real, {0, 2  π}]; 
aa = 
              Random[Real, {0, 
                  2  π}]; 
trajectory3  =  
              Table[{}, {stepsper * cycles}]; 
RungeAll[
              0, a, aa, stepsper, cycles]; 
Do[
              trajectory3[[k]]  =  {Mod[trajectory2[[k]][[1]] - π, 
                    2  π], trajectory2[[k]][[2]]}
, {k, 1, 
                stepsper * cycles}]; 
trajec3 = 
              Join[trajec3, 
                trajectory3];}, {γ}]; 
traject4 = 
          Table[trajec3[[k + killtransients]], {k, 1, 
              cycles*stepsper - killtransients}]; 
ListPlot[
          traject4, PlotStyle -> PointSize[0.002], 
          AxesLabel -> {θ, ω}];
]

Explaination and definitions

This function plots the solution to the set of differential equations defined by parameter values input into this function in terms of θ vs ω.  Note: the two parameters, drivephase and drivefreq, are set equal to, 0 and 2/3, respectively in all plots to reduce the number of inputs required to plot said function, also initial conditions for θ and ω are randomly selected for the same reasons.  
	Local Input Values:
					driverforce - value which sets the amplitude of the driving force function
					dampening -  value which sets the linear dampening constant 
					numofsolutions - value which sets the number of solutions to calculate for this given parameter set
					killtransients - value which sets the number of points to chop from the beginning of the final solution set
	Hidden Local Values:
					stepspercycle -   value which sets the number of steps to take per drivecycle ( where drivecycles are a measure of elapsed time)
					divecycles - value which sets the total number of drivecycles to integrate forward in time 
	Global Parameters Called: 
					there are no direct references to any global parameters previously defined 
	Global Functions Called:	
					SetPendParam - called to set the equations with parameter values input into this function
					RungeAll - called to calculate the values of {θ, ω, ϕ} at each point over the input integration time
					ListPlot - called to output the table of data computed by RungeAll
	Global Parameter Defined:
					a - random real value in the range {0, 2π}which sets the initial value of θ
					aa - random real value in the range {0, 2π}which sets the initial value of ω
					trajec3 - table which contains values for all combined solutions
					trajectory3 - table which contains values for one solution, then written over to contain the second solution, etc.
					traject4 - table which contains all information from trajectory3, except the transients are killed by subtracting the first killtransients points

Examples

As an example, we can call this function with driverforce = 1.5, dampening = 0.5, stepspercycle = 150, cycles = 140, numofsolutions = 1, killtransients = 1400,  and get a listplot of one solution of θ vs ω over 140 drivecycles, dropping the first 1400 terms, for a total of 19600 total points evaluated, and thus we aquire a plot the solution for the parameter set {1.5, 0.5, 2/3, 0}

PlotPhase[1.5, 0.5, 150, 140, 1, 10000]

As a secondary example, we can call this function with driverforce = 1.03, dampening = 0.5, stepspercycle = 100, cycles = 2000, numofsolutions = 1, killtransients = 10000,  and get a listplot of one solution of θ vs ω over 2000 drivecycles, dropping the first 10000 terms, for a total of 190000 total points evaluated, and thus we aquire a plot the solution for the parameter set {1.03, 0.5, 2/3, 0}

PlotPhase[1.03, 0.5, 100, 200, 1, 10000]

Note: As long as the killtransients term is large enough, what is output is an attractor for the system.  The plot represents where any representative point will end up after a long enough integration period.

Poincare Section Functions

Code which defines the function

Poincare[driverforce_, driverfreq_, dampening_, driverphase_, 
        initialt_, initialθ_, initialω_, stepspercycle_, 
        drivecycles_, 
        opts__] := 
Module[{stepsper = stepspercycle, 
          cycles = drivecycles}, 
SetPendParam[driverforce, 
          driverfreq, dampening, 
          driverphase]; 
poincare  =  
          Table[{}, {cycles - ((cycles/10))}]; 
RungeOncePerCycle[initialt, initialθ, 
          initialω, stepsper, cycles]; 
Do[
          poincare[[k]]  =  {Mod[trajectory[[k + 
                          (cycles/10)]][[1]] - π, 
                2  π], trajectory[[k + 
                      (cycles/10)]][[2]]}
, {k, 1, 
            cycles - ((cycles/10))}]; 
ListPlot[
          poincare, PlotStyle -> PointSize[0.0025], 
          opts];
]

Explaination and definitions

This function plots the solution to the set of differential equations defined by parameter values input into this function in terms of θ vs ω.  Note: the two parameters, drivephase and drivefreq, are set equal to, 0 and 2/3, respectively in all plots to reduce the number of inputs required to plot said function, also initial conditions for θ and ω are randomly selected for the same reasons.  
	Local Input Values:
					driverforce - value which sets the amplitude of the driving force function
					drivefreq - value which sets the frequency of the driving force function
					dampening -  value which sets the linear dampening constant 
					driverphase -   value which sets the initial phase of the driving force function
					initialt - value which sets the initial value of t to solve from
					initialθ - value which sets the initial value of θ to solve from
					initialω -  value which sets the initial value of ω to solve from
					opts - options which you want to affect the output listplot from the existing mathematica listplot options
	Hidden Local Values:
					stepspercycle -   value which sets the number of steps to take per drivecycle ( where drivecycles are a measure of elapsed time)
					divecycles - value which sets the total number of drivecycles to integrate forward in time 
	Global Parameters Called: 
					trajectory - is called on from RungeOncePerCycle
	Global Functions Called:	
					SetPendParam - called to set the equations with parameter values input into this function
					RungeOncePerCycle - called to calculate the values of {θ, ω, ϕ} at each point over the input integration time
					ListPlot - called to output the table of data computed by RungeOncePerCycle
	Global Parameter Defined:
					poincare - A table created to hold values of θ and ω, one set per drivecycle, for every cycle integrated over
Note: This function automatically drops the first 10% of the points originally saved, this is done to kill transient solutions.

Examples

As an example, we can call this function with driverforce = 2.3, drivefreq = 1, dampening = 0.05, drivephase = 0, initialt = 0, initialθ =0, initialω = 1, stepspercycle = 100, cycles = 10000,   and get a listplot of one solution of θ vs ω over 10000 drivecycles, dropping the first 1000 terms, for a total of 9000 total points evaluated, and thus we aquire a plot the solution for the parameter set {2.3, 0.05, 1, 0}.  Solutions are also shown for parameter sets {1.5, 0.5, 2/3, 0}with all the same inputs as stated above,  and {1.5, 0.05, 2/3, 0} with a total of 45000 points plotted for the last example, with all other inputs the same.  The third plot shown is just a closeup view of the folding fractal nature of the attractors considered here.

Poincare[2.3, 3/3, 0.05, 0, 0, 1, 1, 100, 10000, 
      PlotRange -> All]

Poincare[1.5, 2/3, 0.5, 0, 1, 1, 100, 10000, 
      PlotRange -> All]

Poincare[1.5, 2/3, 0.5, 0, 1, 1, 100, 20000, 
      PlotRange -> {{4.35, 4.65}, {1.5, 1.75}}]

Poincare[1.5, 2/3, 0.05, 0, 1, 1, 100, 50000, 
      PlotRange -> All]

Creating Animations of Poincare Plots as one parameter is varied over a range

We can make an animation by making successive Poincare Plots of one changing parameter value over some range, and setting them in sequence.  Essentially just Do[] the Poincare function over a range of one of the parameters.

Do[Poincare[1.5, 2/3, 0.25, c, 0, 1, 1, 100, 25000, 
        PlotRange -> All], {c, 0, 2  π, (π/16)}]

Do[Poincare[1.5, 2/3, 0.06, c, 0, 1, 1, 100, 25000, 
        PlotRange -> All], {c, 0, 2  π, (π/16)}]

Bifurcations over Dampening Parameter

Code which defines the function

PendDampBifurcation[initialdamp_, finaldamp_, stepspercycle_, 
        drivecycles_, iterations_, 
        driveforce_] := 
Module[{}, 
θbifurcation =  {{0, 0}}; ωbifurcation  = {{0, 
              0}}; 
Do[
{α = 
              Random[Real, {initialdamp, 
                  finaldamp}]; 
SetPendParam[driveforce, 
              2/3, α, 0]; 
randomθ = 
              Random[Real, {-π, π}]; 
            randomω = 
              Random[Real, {-π, π}]; 
RungeOncePerCycle[0, randomθ, randomω, 
              stepspercycle, 
              drivecycles]; 
βθ = 
              Mod[trajectory[[drivecycles]][[1]] - π, 
                2  π]; βω = trajectory[[drivecycles]][[2]]; 
singleθ  = {{α, βθ}}; 
            singleω  = {{α, βω}}; 
θbifurcation  =  
              Join[singleθ, θbifurcation]; ωbifurcation  =  Join[
                singleω, ωbifurcation];}
, {iterations}]; 
p1 = 
          ListPlot[θbifurcation, DisplayFunction -> Identity]; 
        p2 = ListPlot[ωbifurcation, 
            DisplayFunction -> Identity]; 
Show[
          GraphicsArray[{p1, p2}]];
]

Explaination and definitions

This function plots the bifurcation diagram for a set of solutions to the set of differential equations defined by parameter values and range input into this function, in terms of θ vs ω.
Essentially, it integrates the differential equations given by 3 set parameters, while giving the fourth a range which are randomly selected for each iteration.  It integrates each parameter set to some given time ϕ = drivecycles * 3π.  At this time, the value of θ and ω are saved in lists with the random parameter value used to obtain those solution values.  This is done a given number, iterations, times to produce a table of lists of parameter values and solutions values.  Note: the two parameters, drivephase and drivefreq, are set equal to, 0 and 2/3, respectively in all plots to reduce the number of inputs required to plot said function, also initial conditions for θ and ω are randomly selected for the same reasons.  
	Local Input Values:
					initialdamp -  value which sets the initial linear dampening constant, helping set the overall range
					finaldamp -  value which sets the final linear dampening constant, helping set the overall range 
					stepspercycle -   value which sets the number of steps to take per drivecycle ( where drivecycles are a measure of elapsed time)
					drivecycles - value which sets the total number of drivecycles to integrate forward in time 
					iterations - value which sets the total number of solutions to calculate, and thus how many points to plot overall
					driverforce - value which sets the amplitude of the driving force function
	Hidden Local Values:
					there are no hidden local variables defined
	Global Parameters Called: 
					trajectory - is called on from RungeOncePerCycle
	Global Functions Called:	
					SetPendParam - called to set the equations with parameter values input into this function
					RungeOncePerCycle - called to calculate the values of {θ, ω, ϕ} at each point over the input integration time
					Join - called to connect each set of parameter value used vs solution value at end of integration
					ListPlot - called to output the table of data computed by RungeOncePerCycle, output is hidden by displayfunction-> identity
					Show - called to output a table of the two hidden Listplot commands used
	Global Parameter Defined:
					θbifurcation - A table of lists, each list contains the dampening value used and the solution to θ[t] after integrating over the designated drivecycle 
					ωbifurcation - A table of lists, each list contains the dampening value used and the solution to ω[t] after integrating over the designated drivecycle 
					α - defined to be the value of the randomly selected dampening parameter defining the equation set currently being integrated over
					randomθ - defined to be a randomly selected initial condition for θ
					randomω - defined to be a randomly selected initial condition for ω
					βθ - defined to be the value of θ[t] after integrating over the designated number of drivecycles, each iteration overwrites the previous
					βω - defined to be the value of ω[t] after integrating over the designated number of drivecycles, each iteration overwrites the previous
					singleθ - single elements which build θbifurcation by Join every iteration
					singleω - single elements which build ωbifurcation by Join every iteration
					p1 - defined as the Listplot of θbifurcation, hidden at original calling, displayed as table with p2
					p2 - defined as the Listplot of ωbifurcation, hidden at original calling, displayed as table with p1

Examples

As an example, we can call this function with an initialdamp = 0.05, finaldamp = 0.5, with stepspercycle = 100, drivecycles = 75, iterations 16000, and driverforce = 1.5.  We get a table of listplots of sets of parameter values and solution values of θ and ω over 16000 different parameter sets.  We thus have a plot the possible solution values at ϕ = 75*3π vs the parameter value used to obtain the given solution, this is done for θ and ω giving two plots, over the range {0.05 ≤ b ≤ 0.5}  for the parameter set {1.5, {0.05 ≤ b ≤ 0.5}, 2/3, 0}.  Solutions are also shown for parameter sets {1.0, {0.05 ≤ b ≤ 0.5}, 2/3, 0}with all the same inputs as stated above except half as many iterations.

PendDampBifurcation[0.05, 0.5, 100, 75, 16000, 1.5, 
      PlotStyle -> PointSize[0.002]]

PendDampBifurcation[0.05, 0.5, 100, 75, 8000, 1.0, 
      PlotStyle -> PointSize[0.002]]

Bifurcations over Drive Force Parameter

Code which defines the function

PendForceBifurcation[initialforce_, finalforce_, slices_, cycles_, 
        steps_, dampening_] := 
Module[{}, 
θbifurcation = {{0, 
              0}}; ωbifurcation  =  {{0, 
              0}}; 
Do[
{α = 
              Random[Real, {initialforce, 
                  finalforce}]; 
SetPendParam[α, 
              2/3, dampening, 0]; 
randomθ = 
              Random[Real, {-π, π}]; 
            randomω = 
              Random[Real, {-π, π}]; 
Runge[0, 
              randomθ, randomω, slices, 
              cycles]; 
singleθ  = {{α, βθ}}; 
            singleω  = {{α, βω}}; 
βθ = 
              Mod[trajectory[[cycles]][[1]] - π, 
                2  π]; βω = trajectory[[cycles]][[2]]; 
θbifurcation  =  Join[singleθ, θbifurcation]; ωbifurcation  =  
              Join[singleω, ωbifurcation];}
, {steps}]; 
p1 = 
          ListPlot[θbifurcation, DisplayFunction -> Identity]; 
        p2 = ListPlot[ωbifurcation, 
            DisplayFunction -> Identity]; 
Show[
          GraphicsArray[{p1, p2}]];
]

Explaination and definitions

This function plots the bifurcation diagram for a set of solutions to the set of differential equations defined by parameter values and range input into this function, in terms of θ vs ω.
Essentially, it integrates the differential equations given by 3 set parameters, while giving the fourth a range which are randomly selected for each iteration.  It integrates each parameter set to some given time ϕ = drivecycles * 3π.  At this time, the value of θ and ω are saved in lists with the random parameter value used to obtain those solution values.  This is done a given number, iterations, times to produce a table of lists of parameter values and solutions values.  Note: the two parameters, drivephase and drivefreq, are set equal to, 0 and 2/3, respectively in all plots to reduce the number of inputs required to plot said function, also initial conditions for θ and ω are randomly selected for the same reasons.  
	Local Input Values:
					initialforce -  value which sets the initial linear dampening constant, helping set the overall range
					finalforce -  value which sets the final linear dampening constant, helping set the overall range 
					stepspercycle -   value which sets the number of steps to take per drivecycle ( where drivecycles are a measure of elapsed time)
					drivecycles - value which sets the total number of drivecycles to integrate forward in time 
					iterations - value which sets the total number of solutions to calculate, and thus how many points to plot overall
					dampening - value which sets the amplitude of the driving force function
	Hidden Local Values:
					there are no hidden local variables defined
	Global Parameters Called: 
					trajectory - is called on from RungeOncePerCycle
	Global Functions Called:	
					SetPendParam - called to set the equations with parameter values input into this function
					RungeOncePerCycle - called to calculate the values of {θ, ω, ϕ} at each point over the input integration time
					Join - called to connect each set of parameter value used vs solution value at end of integration
					ListPlot - called to output the table of data computed by RungeOncePerCycle, output is hidden by displayfunction-> identity
					Show - called to output a table of the two hidden Listplot commands used
	Global Parameter Defined:
					θbifurcation - A table of lists, each list contains the dampening value used and the solution to θ[t] after integrating over the designated drivecycle 
					ωbifurcation - A table of lists, each list contains the dampening value used and the solution to ω[t] after integrating over the designated drivecycle 
					α - defined to be the value of the randomly selected dampening parameter defining the equation set currently being integrated over
					randomθ - defined to be a randomly selected initial condition for θ
					randomω - defined to be a randomly selected initial condition for ω
					βθ - defined to be the value of θ[t] after integrating over the designated number of drivecycles, each iteration overwrites the previous
					βω - defined to be the value of ω[t] after integrating over the designated number of drivecycles, each iteration overwrites the previous
					singleθ - single elements which build θbifurcation by Join every iteration
					singleω - single elements which build ωbifurcation by Join every iteration
					p1 - defined as the Listplot of θbifurcation, hidden at original calling, displayed as table with p2
					p2 - defined as the Listplot of ωbifurcation, hidden at original calling, displayed as table with p1

Examples

PendForceBifurcation[0.9, 1.75, 100, 75, 8000, 0.5, 
      PlotStyle -> PointSize[0.002]]

PendForceBifurcation[0.9, 1.75, 100, 75, 8000, 0.2, 
      PlotStyle -> PointSize[0.002]]

Calculating Lyapunov Exponents Via Comparing Neighboring Trajectories

Calculating Solutions For θ[t] and Saving Data

Code which defines the function

CalculateTheta[driverforce_, dampening_, stepspercycle_, drivecycles_, 
        initialθ_, 
        initialω_] := 
Module[{stepsper = 
            stepspercycle, 
          cycles = drivecycles}, 
SetPendParam[driverforce, 
          2/3, dampening, 0]; 
trajectory4  =  
          Table[{}, {stepsper * cycles}]; 
RungeAll[0, 
          initialθ, initialω, stepsper, 
          cycles]; 
Do[
          trajectory4[[k]]  =  trajectory2[[k]][[1]], {k, 
            1, stepsper * cycles}]; 
Return[
          trajectory4];
]

Explaination and definitions

This function calculates the solution to the set of differential equations defined by parameter values input into this function in terms of θ vs t.  Note: the two parameters, drivephase and drivefreq, are set equal to, 0 and 2/3, respectively in all plots to reduce the number of inputs required to plot said function, also initial conditions for θ and ω are randomly selected for the same reasons.  
	Local Input Values:
					driverforce - value which sets the amplitude of the driving force function
					dampening -  value which sets the linear dampening constant 
					initialθ - value which sets the initial value of θ to solve from
					initialω -  value which sets the initial value of ω to solve from
	Hidden Local Values:
					stepspercycle -   value which sets the number of steps to take per drivecycle ( where drivecycles are a measure of elapsed time)
					divecycles - value which sets the total number of drivecycles to integrate forward in time 
	Global Parameters Called: 
					trajectory2 is called on from RungeAll
	Global Functions Called:	
					SetPendParam - called to set the equations with parameter values input into this function
					RungeAll - called to calculate the values of {θ, ω, ϕ} at each point over the input integration time
	Global Parameter Defined:
					trajectory4 - table which contains values for the solution, contains entries for every step over each drivecycle, is returned when function gets called

Examples

As an example we set driveforce = 1.5, dampening = 0.5, stepspercycle = 100, drivecycles = 10, and initialθ = initialω = 1.  When we call this function with those inputs we are output a list of data containing the value of θ and the value of t over the given range.  Instead of showing the list of data, ive shown the listplot of the data to make it more clear.

CalculateTheta[1.5, 0.5, 100, 10, 1, 1];
ListPlot[%];

Calculating Solutions For ω[t] and Saving Data

Code which defines the function

CalculateOmega[driverforce_, dampening_, stepspercycle_, drivecycles_, 
        initialθ_, 
        initialω_] := 
Module[{stepsper = 
            stepspercycle, 
          cycles = drivecycles}, 
SetPendParam[driverforce, 
          2/3, dampening, 0]; 
trajectory5  =  
          Table[{}, {stepsper * cycles}]; 
RungeAll[0, 
          initialθ, initialω, stepsper, 
          cycles]; 
Do[
          trajectory5[[k]]  =  trajectory2[[k]][[2]], {k, 
            1, stepsper * cycles}]; 
Return[
          trajectory5];
]

Explaination and definitions

This function calculates the solution to the set of differential equations defined by parameter values input into this function in terms of ω vs t.  Note: the two parameters, drivephase and drivefreq, are set equal to, 0 and 2/3, respectively in all plots to reduce the number of inputs required to plot said function, also initial conditions for θ and ω are randomly selected for the same reasons.  
	Local Input Values:
					driverforce - value which sets the amplitude of the driving force function
					dampening -  value which sets the linear dampening constant 
					initialθ - value which sets the initial value of θ to solve from
					initialω -  value which sets the initial value of ω to solve from
	Hidden Local Values:
					stepspercycle -   value which sets the number of steps to take per drivecycle ( where drivecycles are a measure of elapsed time)
					divecycles - value which sets the total number of drivecycles to integrate forward in time 
	Global Parameters Called: 
					trajectory2 is called on from RungeAll
	Global Functions Called:	
					SetPendParam - called to set the equations with parameter values input into this function
					RungeAll - called to calculate the values of {θ, ω, ϕ} at each point over the input integration time
	Global Parameter Defined:
					trajectory5 - table which contains values for the solution, contains entries for every step over each drivecycle, is returned when function gets called

Examples

As an example we set driveforce = 1.5, dampening = 0.5, stepspercycle = 100, drivecycles = 10, and initialθ = initialω = 1.  When we call this function with those inputs we are output a list of data containing the value of θ and the value of t over the given range.  Instead of showing the list of data, ive shown the listplot of the data to make it more clear.

CalculateOmega[1.5, 0.5, 100, 10, 1, 1];
ListPlot[%];

Comparing Solutions For θ[t] With Nearby Initial Conditions

Code which defines the function

CompareTheta[driverforce_, dampening_, stepspercycle_, drivecycles_, 
        initialθ_, initialω_, 
        delta_] := 
Module[{stepsper = stepspercycle, 
          cycles = drivecycles, δ = 
            delta}, 
solution1  =  
          CalculateTheta[driverforce, dampening, stepsper, cycles, 
            initialθ, 
            initialω]; 
solution2  =  
          CalculateTheta[driverforce, dampening, stepsper, cycles, 
            initialθ + δ, 
            initialω + δ]; 
difference  =  
          Table[{0}, { .9*stepsper*cycles}]; 
Do[
          difference[[i]]  =  
            Log[Abs[solution2[[i +  .1*stepsper*cycles]] - 
                  solution1[[i +  .1*stepsper*cycles]]]], {i, 
            1,  .9*stepsper*cycles}]; 
differenceplot  = 
          ListPlot[difference, 
            DisplayFunction -> Identity];
]

Explaination and definitions

This function calculates the solutions to the set of differential equations defined by parameter values input into this function in terms of θ vs t, but also then does the same for the same parameter set with initial conditions varied slightly.  The difference between these two trajectories is then calculated over the entire range.  A listplot is the output, which is hidden.  Note: the two parameters, drivephase and drivefreq, are set equal to, 0 and 2/3, respectively in all plots to reduce the number of inputs required to plot said function, also initial conditions for θ and ω are randomly selected for the same reasons.  
	Local Input Values:
					driverforce - value which sets the amplitude of the driving force function
					dampening -  value which sets the linear dampening constant 
					initialθ - value which sets the initial value of θ to solve from
					initialω -  value which sets the initial value of ω to solve from
					delta - value which sets the initial offset of initial conditions in the two solutions
	Hidden Local Values:
					stepspercycle -   value which sets the number of steps to take per drivecycle ( where drivecycles are a measure of elapsed time)
					divecycles - value which sets the total number of drivecycles to integrate forward in time 
	Global Parameters Called: 
					no global parameters are called on explicitly in this function
	Global Functions Called:	
					CalculateTheta is called on twice to produce solutions to the set of equations
	Global Parameter Defined:
					solution1 - is defined to be the solution to the set of differential equations, in terms of θ and t.
					solution2 - is defined to be the solution to the set of differential equations, with offset initial conditions, in terms of θ and t.
					difference - is defined to be the difference in the two solutions at ever point integrated
					differenceplot - is defined to be the hidden listplot output of this function.

Examples

As an example, we can call this function with driveforce = 1.24, dampening 0.5, stepspercycle = 100, drivecycles = 10, initialθ = initialω = 1, delta = 0.01

CompareTheta[1.24, 0.5, 100, 1, 1, 
        1,  .01];
Show[differenceplot, 
        DisplayFunction -> $DisplayFunction];

Comparing Solutions For ω[t] With Nearby Initial Conditions

Code which defines the function

CompareOmega[driverforce_, dampening_, stepspercycle_, drivecycles_, 
        initialθ_, initialω_, 
        delta_] := 
Module[{stepsper = stepspercycle, 
          cycles = drivecycles, δ = 
            delta}, 
solution3  =  
          CalculateOmega[driverforce, dampening, stepsper, cycles, 
            initialθ, 
            initialω]; 
solution4  =  
          CalculateOmega[driverforce, dampening, stepsper, cycles, 
            initialθ + δ, 
            initialω + δ]; 
difference2  =  
          Table[{0}, { .9*stepsper*cycles}]; 
Do[
          difference2[[i]]  =  
            Log[Abs[solution4[[i +  .1*stepsper*cycles]] - 
                  solution3[[i +  .1*stepsper*cycles]]]], {i, 
            1,  .9*stepsper*cycles}]; 
differenceplot2  = 
          ListPlot[difference2, 
            DisplayFunction -> Identity];
]

Explaination and definitions

This function calculates the solutions to the set of differential equations defined by parameter values input into this function in terms of ω vs t, but also then does the same for the same parameter set with initial conditions varied slightly.  The difference between these two trajectories is then calculated over the entire range.  A hidden listplot is the output.  Note: the two parameters, drivephase and drivefreq, are set equal to, 0 and 2/3, respectively in all plots to reduce the number of inputs required to plot said function, also initial conditions for θ and ω are randomly selected for the same reasons.  
	Local Input Values:
					driverforce - value which sets the amplitude of the driving force function
					dampening -  value which sets the linear dampening constant 
					initialθ - value which sets the initial value of θ to solve from
					initialω -  value which sets the initial value of ω to solve from
					delta - value which sets the initial offset of initial conditions in the two solutions
	Hidden Local Values:
					stepspercycle -   value which sets the number of steps to take per drivecycle ( where drivecycles are a measure of elapsed time)
					divecycles - value which sets the total number of drivecycles to integrate forward in time 
	Global Parameters Called: 
					no global parameters are called on explicitly in this function
	Global Functions Called:	
					CalculateTheta is called on twice to produce solutions to the set of equations
	Global Parameter Defined:
					solution1 - is defined to be the solution to the set of differential equations, in terms of ω and t.
					solution2 - is defined to be the solution to the set of differential equations, with offset initial conditions, in terms of ω and t.
					difference - is defined to be the difference in the two solutions at ever point integrated
					differenceplot - is defined to be the hidden listplot output of this function.

Examples

As an example, we can call this function with driveforce = 1.21, dampening 0.5, stepspercycle = 100, drivecycles = 10, initialθ = initialω = 1, delta = 0.01.  We can show the hidden listplot with the show command

CompareOmega[1.21, 0.5, 100, 1, 1, 
        1,  .01];
Show[differenceplot2, 
        DisplayFunction -> $DisplayFunction];

Calculating the Exponential Rate of Divergence of Initially Close Trajectories

Code which defines the function

Lyapunov[driverforce_, dampening_, stepspercycle_, drivecycles_,  
        initialθ_, initialω_,  
        delta_] := 
Module[{stepsper = stepspercycle, 
          cycles = drivecycles, δ = 
            delta}, 
lyapunov1  =  
          0; 
lyapunov2  =  
          0; 
phasediff1  =  
          Table[0, {cycles}]; 
phasediff2  =  
          Table[0, {cycles}]; 
CompareTheta[driverforce, 
          dampening, stepsper, cycles, initialθ, 
          initialω, δ]; 
CompareOmega[
          driverforce, dampening, stepsper, cycles, initialθ, 
          initialω, δ]; 
Do[
          phasediff1[[j]]  =  difference[[j]], {j, 1, 
            cycles}]; 
Do[
          phasediff2[[j]]  =  difference2[[j]], {j, 1, 
            cycles}]; 
lyapunov1  =  
          D[Fit[phasediff1, {1, x}, x], 
            x]; 
lyapunov2  =  
          D[Fit[phasediff2, {1, x}, x], x];
]

Explaination and definitions

This function calculates the exponential rate of divergence of two solutions to the set of differential equations defined by parameter values input into this function.  The differences between two trajectories is then calculated over the entire range using with CompareTheta, and CompareOmega, we then apply a linear curve fir to this set of data, the slope to this fit is then the rate of divergence of two nearby solutions.  One is calculated for θ and one for ω.  Note: the two parameters, drivephase and drivefreq, are set equal to, 0 and 2/3, respectively in all calculations to reduce the number of inputs required to plot said function.  
	Local Input Values:
					driverforce - value which sets the amplitude of the driving force function
					dampening -  value which sets the linear dampening constant 
					initialθ - value which sets the initial value of θ to solve from
					initialω -  value which sets the initial value of ω to solve from
					delta - value which sets the initial offset of initial conditions in the two solutions
	Hidden Local Values:
					stepspercycle -   value which sets the number of steps to take per drivecycle ( where drivecycles are a measure of elapsed time)
					divecycles - value which sets the total number of drivecycles to integrate forward in time 
	Global Parameters Called: 
					difference1 - is called from the function comparetheta
					difference2 - is called from the function compareomega
	Global Functions Called:	
					CompareTheta is called on to produce solutions to the set of equations
					CompareOmega is called on to produce solutions to the set of equations
	Global Parameter Defined:
					phasediff1 - is defined to be the table containing the differences in the θ direction over the whole range
					phasediff2 - is defined to be the table containing the differences in the ω direction over the whole range
					lyapunov1 - is defined to be the rate of divergence in the θ direction
					lyapunov2 - is defined to be the rate of divergence in the ω direction

Examples

As an example, we can calculate the rate of divergence of the parameter set with driveforce = 0.9, dampening = 0.5, with stepspercycle = 100, and over 1 drivecycle, with initialθ = initialω = 1, and delta = 0.01.  This gives two exponent values, both negative, assuring we have periodic behavior in this region.

Lyapunov[0.9, 0.5, 100, 1, 1, 1,  .01]
{lyapunov1, lyapunov2}

Plot Range of Exponential Rates for a Range of Drive Force Values

Code which defines the function

PlotLyapunov[dampening_, stepspercycle_, drivecycles_,   delta_, 
        totalsamplesize_] := 
Module[{stepsper = 
            stepspercycle, 
          cycles = drivecycles, δ = 
            delta}, 
lyapunovset  =  {{0, 
              0}}; 
lyapunovset2  =  {{0, 
              0}}; 
Do[{randomf = 
              Random[Real, {0.9, 1.8}]; 
randomθ = 
              Random[Real, {-π, π}]; 
            randomω = 
              Random[Real, {-π, π}]; 
Lyapunov[
              randomf, dampening, stepsper, cycles, randomθ, 
              randomω, δ]; 
β1  =  
              lyapunov1; 
β2  =  
              lyapunov2; 
singlepoint  =  {{randomf, β1}}; 
singlepoint2  =  {{randomf, β2}}; 
lyapunovset  =  
              Join[singlepoint, 
                lyapunovset]; 
lyapunovset2  =  
              Join[singlepoint2, 
                lyapunovset2];
}, {totalsamplesize}]; 
lyapunovset  =  
          Drop[lyapunovset,  -1]; 
lyapunovset2  =  
          Drop[lyapunovset2,  -1]; 
p1  =  
          ListPlot[lyapunovset /. Indeterminate -> 0, 
            DisplayFunction -> Identity]; 
p2  =  
          ListPlot[lyapunovset2 /. Indeterminate -> 0, 
            DisplayFunction -> Identity]; 
Show[
          GraphicsArray[{p1, p2}]];
]

Explaination and definitions

This function calculates the exponential rate of divergence of two solutions to the set of differential equations defined by some range of parameter values input into this function.  The differences between two trajectories is then calculated over the entire range using with CompareTheta, and CompareOmega, we then call Lyapunov to calculate the linear slope of the diverging trajectories.  One is calculated for θ and one for ω.  Note: the two parameters, drivephase and drivefreq, are set equal to, 0 and 2/3, respectively in all calculations to reduce the number of inputs required to plot said function, also initialθ and initialω are randomly selected.  
	Local Input Values:
					dampening -  value which sets the linear dampening constant 
					delta - value which sets the initial offset of initial conditions in the two solutions
					totalsamplesize - total number of points to calculate
	Hidden Local Values:
					stepspercycle -   value which sets the number of steps to take per drivecycle ( where drivecycles are a measure of elapsed time)
					divecycles - value which sets the total number of drivecycles to integrate forward in time 
	Global Parameters Called: 
					lyapunov1 - is called from the function Lyapunov
					lyapunov2 - is called from the function Lyapunov
	Global Functions Called:	
					Lyapunov - is called to calculate the rate of divergence of trajectories
	Global Parameter Defined:
					lyapunovset - set which contains the θ lyapunov exponents over the whole range, each value is paired with the drive force value which produced that exponent
					lyapunovset2 - set which contains the θ lyapunov exponents over the whole range, each value is paired with the drive force value which produced that exponent
					randomf - randomly selected drive force value which determines the parameter set to be calculated
					randomθ - defined to be a randomly selected initial condition for θ
					randomω - defined to be a randomly selected initial condition for ω
					singlepoint - set which contains the exponent value as well as the driveforce value which was used, overwritten each iteration
					singlepoint2 - set which contains the exponent value as well as the driveforce value which was used, overwritten each iteration
					β1 - value which contains each θ exponent value, overwritten by next iterations exponent value, etc
					β2 - value which contains each ω exponent value, overwritten by next iterations exponent value, etc
					p1 - name given to listplot of the range of θ Lyapunov exponents 
					p2 - name given to listplot of the range of ω Lyapunov exponents

Examples

As an example, we can calculate the rate of divergence of the parameter set with  dampening = 0.5, with stepspercycle = 100, and over 1 drivecycle, with delta = 0.1, and 500 total points plot.

PlotLyapunov[0.5, 100, 1, 0.1, 500]

Plot Basins of Exponential Divergence Rate Over Range Of Initial Conditions

Code which defines the function

BasinValue[
        slicewidth_, θstart_, θfinish_, ωstart_, ωfinish_, delta_, stepspercycle_, driverforce_, 
        opts__] := 
Module[{}, 
basinθwidth  =  
          Floor[N[((θfinish - θstart)/slicewidth)]]; 
basinωwidth  =  
          Floor[N[((ωfinish - ωstart)/slicewidth)]]; 
totalbasin  =  {Table[
              0, {basinωwidth}]}; 
Do[{oneslice  =  Table[{0}, {basinθwidth}]; 
lyapanov11  =  {{0, 
                  0}}; 
lyapanov22  =  {{0, 
                  0}}; 
Do[{
Lyapunov[
                  driverforce, 0.5, stepspercycle, 1, 
                  delta, γ, η]; 
lyap1  =  {{γ, Re[lyapunov1]}}; 
lyapanov11  =  
                  Join[lyap1, 
                    lyapanov11]; 
lyap2  =  {{γ, 
                      Re[lyapunov2]}}; 
lyapanov22  =  
                  Join[lyap2, 
                    lyapanov22];
}, {γ, θstart, θfinish, 
                slicewidth}]; 
lyapanov11  =  
              Drop[lyapanov11, -1]; 
lyapanov22  =  
              Drop[lyapanov22, -1]; 
Do[{oneslice[[j]] = lyapanov11[[j]][[2]]}, {j, 1, 
                basinθwidth}]; 
totalbasin  =  
              Join[{oneslice}, 
                totalbasin];}, {η, ωstart, ωfinish, 
            slicewidth}]; 
totalbasin  =  
          Drop[totalbasin, -1]; 
ListContourPlot[
          totalbasin, ContourLines -> False, 
          MeshRange -> {{θstart, θfinish}, {ωstart, ωfinish}}, opts];
]

Explaination and definitions

This function calculates the exponential rate of divergence of two solutions to the set of differential equations defined by some parameter values input into this function, over a range of initial conditions.  The differences between two trajectories is then calculated over the entire range using with CompareTheta, and CompareOmega, we then call Lyapunov to calculate the linear slope of the diverging trajectories.  One is calculated for θ and one for ω.  Note: the two parameters, drivephase and drivefreq, are set equal to, 0 and 2/3, respectively in all calculations to reduce the number of inputs required to plot said function, also initialθ and initialω are randomly selected.  
	Local Input Values:
					slicewidth - how often to sample the initial conditions
					θstart - first θ value to input
					θfinish - last θ value to input
					ωstart - first ω value to input
					ωfinish - last ω value to input
					delta - value which sets the initial offset of initial conditions in the two solutions
					driverforce - value of the drive force to input as the parameter set
	Hidden Local Values:
					stepspercycle -   value which sets the number of steps to take per drivecycle ( where drivecycles are a measure of elapsed time)
	Global Parameters Called: 
					lyapunov1 - is called from the function Lyapunov
					lyapunov2 - is called from the function Lyapunov
	Global Functions Called:	
					Lyapunov - is called to calculate the rate of divergence of trajectories
					ListContourPlot - is called to output the results.
	Global Parameter Defined:
					lyapunov11 - set which contains the θ lyapunov exponents over the whole range, each value is paired with the drive force value which produced that exponent
					lyapunov22 - set which contains the θ lyapunov exponents over the whole range, each value is paired with the drive force value which produced that exponent
					lyap1 - set which contains the exponent value as well as the driveforce value which was used, overwritten each iteration
					lyap2 - set which contains the exponent value as well as the driveforce value which was used, overwritten each iteration
					oneslice - one full range of the θ direction, while still just one value for ω
					totalbasin - combination of each oneslice

Examples

BasinValue[ .1, -π, π, -π, π,  .1, 100, 0.9, 
      DisplayFunction -> $DisplayFunction]

Plot Basins of Exponential Divergence Rate Over Range Of Initial Conditions 2

Code which defines the function

BasinCutOff[
        slicewidth_, θstart_, θfinish_, ωstart_, ωfinish_, phaseδ_, stepspercycle_, driverforce_, 
        opts__] := 
Module[{}, 
basinθwidth  =  
          Floor[N[((θfinish - θstart)/slicewidth)]]; 
basinωwidth  =  
          Floor[N[((ωfinish - ωstart)/slicewidth)]]; 
totalbasin  =  {Table[
              0, {basinωwidth}]}; 
Do[{oneslice  =  Table[{0}, {basinθwidth}]; 
lyapanov11  =  {{0, 
                  0}}; 
lyapanov22  =  {{0, 
                  0}}; 
Do[{
Lyapunov[
                  driverforce, 0.5, stepspercycle, 25, 
                  phaseδ, γ, η]; 
lyap1  =  {{γ, Re[lyapunov1]}}; 
lyapanov11  =  
                  Join[lyap1, 
                    lyapanov11]; 
lyap2  =  {{γ, 
                      Re[lyapunov2]}}; 
lyapanov22  =  
                  Join[lyap2, 
                    lyapanov22];
}, {γ, θstart, θfinish, 
                slicewidth}]; 
lyapanov11  =  
              Drop[lyapanov11, -1]; 
lyapanov22  =  
              Drop[lyapanov22, -1]; 
templyap1  = lyapanov11[[η]][[2]]; 
templyap2  = lyapanov22[[η]][[2]]; 
If[
              templyap1 > templyap2,  nw  =  templyap1, 
              nw  =  templyap2]; 
If[
              templyap1 < templyap2,  nw2  =  templyap2, 
              nw2  =  
                templyap1]; 
Do[{oneslice[[j]] = 
                  If[nw > 0, 1, 0]}, {j, 1, 
                basinθwidth}]; 
totalbasin  =  
              Join[{oneslice}, 
                totalbasin];}, {η, ωstart, ωfinish, 
            slicewidth}]; 
totalbasin  =  
          Drop[totalbasin, -1]; 
ListContourPlot[
          totalbasin, ContourLines -> False, 
          
          MeshRange -> {{θstart, θfinish}, {ωstart, ωfinish}}, opts];
]

Explaination and definitions

This function calculates the exponential rate of divergence of two solutions to the set of differential equations defined by some parameter values input into this function, over a range of initial conditions.  The differences between two trajectories is then calculated over the entire range using with CompareTheta, and CompareOmega, we then call Lyapunov to calculate the linear slope of the diverging trajectories.  One is calculated for θ and one for ω.  Note: the two parameters, drivephase and drivefreq, are set equal to, 0 and 2/3, respectively in all calculations to reduce the number of inputs required to plot said function, also initialθ and initialω are randomly selected.  
	Local Input Values:
					slicewidth - how often to sample the initial conditions
					θstart - first θ value to input
					θfinish - last θ value to input
					ωstart - first ω value to input
					ωfinish - last ω value to input
					delta - value which sets the initial offset of initial conditions in the two solutions
					driverforce - value of the drive force to input as the parameter set
	Hidden Local Values:
					stepspercycle -   value which sets the number of steps to take per drivecycle ( where drivecycles are a measure of elapsed time)
	Global Parameters Called: 
					lyapunov1 - is called from the function Lyapunov
					lyapunov2 - is called from the function Lyapunov
	Global Functions Called:	
					Lyapunov - is called to calculate the rate of divergence of trajectories
					ListContourPlot - is called to output the results.
	Global Parameter Defined:
					lyapunov11 - set which contains the θ lyapunov exponents over the whole range, each value is paired with the drive force value which produced that exponent
					lyapunov22 - set which contains the θ lyapunov exponents over the whole range, each value is paired with the drive force value which produced that exponent
					lyap1 - set which contains the exponent value as well as the driveforce value which was used, overwritten each iteration
					lyap2 - set which contains the exponent value as well as the driveforce value which was used, overwritten each iteration
					oneslice - one full range of the θ direction, while still just one value for ω
					totalbasin - combination of each oneslice

Examples

BasinCutOff[ .1, -2 π, 2  π, -2 π, 2  π,  .1, 
      100, 1.5, DisplayFunction -> $DisplayFunction]