Shaken Washboard, or the Driven Damped Pendulum.
Version 2.1.  New features since version 1 are
  * Parameter nT, for number of Transients .This many periods are integrated numerically but not plotted.
  * A dot in the theta, dtheta/dt plane is plotted after each driving each period.
  * Only the Poincare map, in the theta, dtheta/dt plane is plotted. For this plot, theta is put in the interval [-Pi, pi).
  
This is a modification of a Mathematica notebook written by my son, Michael Swift, for a junior level mechanics course in the physics department at Rice University.  
This is also a mathematical model of the Actual Shaken Washboard (ASW), built by Mechanical Engineering students Doug LaMaster, Jeremy Petak, and Sam Zerbib, who took this class in Fall 2012.
That machine has 
L = 3 (inches)
H = 1/8 (inch)
A = 0.9 (inches)
Those parameters give the period of small oscillations to be almost exactly 1 second.
The parameter a (called A-bar in Homework 2 from Fall 2023) is proportional to A/L, and is fixed at about 0.67 for the ASW.
The motor can vary the driving period between about 0.6 sec and 1.6 sec.  The parameter ω is (1 second)/(driving period).  

To run this program, the user defines the dimensionless parameters
c, a, and ω, as well as the number of driving periods plotted, then runs the cell.

Define a function that solves the system, with the parameters c, ω0, etc. to be defined later.

Clear[θ]
xInIntervalab[x_, a_, b_]:= Module[{s},
s = Floor[(x-a)/(b-a)]*1.; (* the * 1. makes this a floating point number *)
x - s(b-a)];

runPendulum2 := Module[{tFinal = 2 Pi nP/ω, tTransient = 2Pi nT/ω},
Print["Parameters are {a, c, ω} = ", {a, c, ω}];
Print["Initial conditions are {θ, θ'} = ", {θ0, v0}];
sol=NDSolve[{θ''[t]+c *θ'[t]+Sin[θ[t]]==a*ω^2*Cos[ω*t],θ[0]== θ0,θ'[0]==v0},θ,{t,0,tTransient}];
θfinal = θ[tTransient]/.sol[[1]];
θ0 = xInIntervalab[θfinal, -Pi, Pi]; (* new initial conditions after the transients *)
v0 = θ'[tTransient]/.sol[[1]]; (* new initial conditions *)
Print["The new initial conditions are the final state after nT = ", nT, " driving periods:"];
Print["New Initial {θ, θ'} is ", {θ0, v0}];
sol=NDSolve[{θ''[t]+c *θ'[t]+Sin[θ[t]]==a*ω^2*Cos[ω*t],θ[0]== θ0,θ'[0]==v0},θ,{t,0,tFinal}];
(*grid2Pi[θMin_, θMax_] := (* puts grid lines at multiples of 2 Pi *)
Table[2 Pi i, {i, Floor[θMin/(2Pi)] , Ceiling[θMax/(2 Pi)] }];*)
θfinal = θ[tFinal]/.sol[[1]];
θfinalInInterval = xInIntervalab[θfinal, -Pi, Pi];
vfinal = θ'[tFinal]/.sol[[1]];
Print["Final {θ, θ'} is ", {θfinal, vfinal}];
If[θfinal>= Pi || θfinal < -Pi,
Print["Final {θ, θ'} is ", {θfinalInInterval, vfinal}, ", when θ is put in the interval [-π, π)"]
];
poincareList = Table[{θ[tt], θ'[tt]}/.sol[[1]], {tt, 0, tFinal, 2Pi/ω}];
p1 = Plot[θ[tt]/.sol[[1]], {tt,0,tFinal}, AxesLabel ->{"ωt/(2π)", "θ"}, Ticks ->{Table[{2Pi n/ω, n}, {n, 0, nP}], Table[n Pi, {n, -10,10}]}];
p2 = ParametricPlot[{θ[tt], θ'[tt]}/.sol[[1]],{tt,0,tFinal}, AspectRatio -> 1, AxesLabel ->{"θ", "dθ/dt"}, Epilog -> Point[poincareList]];
p3 = ListPlot[xInIntervalab[poincareList, -Pi, Pi], Frame -> True, FrameLabel->{"θ", "dθ/dt"}, PlotRange -> All];
GraphicsRow[{p1, p2, p3}, ImageSize -> 800]
]

Now define the parameters, and run the program.  Note that the value of a = 0.67 is not as interesting as larger a.  We will keep ω fixed at 0.7 and increase A-bar (= a) to see the transition to chaos.

a=5 Pi/7 * 0.9/3;  (* A-bar for the Actual Shaken Washboard, approximately 0.67*)
a = 1;  
ω=.7;  (* Driving frequency devided by natural frequency.  About 0.6 to 1.5 for the ASW *)
c = .2;
nT = 200; (* number of Transients: the number of driving periods integrated before plotting *)
nP = 3; (* number of driving periods Plotted *)
{θ0, v0} = {1,1.5};
runPendulum2

poincareList

Here is a copy of that previous cell.
Now you can use that final state as the initial condition for another run and approximate the periodic solution.

Now change a to 1.3.  Note that this is an “asymmetric solution” that goes higher on the left.

a=5 Pi/7 * 0.9/3;  (* A-bar for the Actual Shaken Washboard, approximately 0.67*)
a = 1.3; 
ω=.7;  (* Driving frequency devided by natural frequency.  About 0.6 to 1.5 for the ASW *)
c = .2;
nT = 20; (* number of Transients: the number of driving periods integrated before plotting *)
nP = 2; (* number of driving periods plotted *)
{θ0, v0} = {-2,1};
runPendulum2

Here is a trick to see the other, symmetrically related Initial condition.  Let nT to be a “half integer” and then change the IC to the negative of the initial condition it uses, and change nT back to an integer (possibly 0)

a=5 Pi/7 * 0.9/3;  (* A-bar for the Actual Shaken Washboard, approximately 0.67*)
a = 1.3; 
ω=.7;  (* Driving frequency devided by natural frequency.  About 0.6 to 1.5 for the ASW *)
c = .2;
nT = 20.5; (* number of Transients: the number of driving periods integrated before plotting *)
nP = 2; (* number of driving periods plotted *)
{θ0, v0} = {-2,1};
runPendulum2

a=5 Pi/7 * 0.9/3;  (* A-bar for the Actual Shaken Washboard, approximately 0.67*)
a = 1.3; 
ω=.7;  (* Driving frequency devided by natural frequency.  About 0.6 to 1.5 for the ASW *)
c = .2;
nT = 0; (* number of Transients: the number of driving periods integrated before plotting *)
nP = 2; (* number of driving periods plotted *)
{θ0, v0} = -{2.462895311116626`,-0.8484923303021182`}; (* the negative of the previous new IC *)
runPendulum2

Continue in this manner with a = 1.39 to see a period doubling.
Then try a = 1.394 and a = 1.396  (It is hard to tell if the attractor of the last one is period 4 or Period 8.)

a=5 Pi/7 * 0.9/3;  (* A-bar for the Actual Shaken Washboard, approximately 0.67*)
a = 1.394; 
ω=.7;  (* Driving frequency devided by natural frequency.  About 0.6 to 1.5 for the ASW *)
c = .2;
nT = 400; (* number of Transients: the number of driving periods integrated before plotting *)
nP = 4; (* number of driving periods plotted *)
{θ0, v0} = {-2,1.15};
runPendulum2

Now go back to a = 1.4, and play with different ICs (Initial Conditions).  This system has multiple attractors, depending on the initial conditions.

a=5 Pi/7 * 0.9/3;  (* A-bar for the Actual Shaken Washboard, approximately 0.67*)
a = 1.4; 
ω=.7;  (* Driving frequency devided by natural frequency.  About 0.6 to 1.5 for the ASW *)
c = .2;
nT = 0; (* number of Transients: the number of driving periods integrated before plotting *)
nP = 6; (* number of driving periods plotted *)
{θ0, v0} = {0, 1}; (* try {-2,1} and {0,1} *)
runPendulum2

a=5 Pi/7 * 0.9/3;  (* A-bar for the Actual Shaken Washboard, approximately 0.67*)
a =2;
ω=.7;  (* Driving frequency devided by natural frequency.  About 0.6 to 1.5 for the ASW *)
c = .1;
nT = 100; (* number of Transients: the number of driving periods integrated before plotting *)
nP = 100; (* number of driving periods plotted *)
{θ0, v0} = {0, 1}; (* try {-2,1} and {0,1} *)
runPendulum2

Here is a simulation of the experiment shown in the youtube video.  This has nT = 0, and initially nP = 20.
But see what happens when you make  nP = 25.  (The nP = 20 and nP25 plots are shown below this, for comparison.)
The notebook ShakenWashboard3.nb does a separate numerical integration for each driving period.  That way, you can make nT and nP large without having conflicting numerical results.  You can never really trust numerical integration of the ODE, but it is not good to make the time interval too large, since the algorithm tries to control the accuracy at the final time.

a=5 Pi/7 * 0.9/3;  (* A-bar for the Actual Shaken Washboard, approximately 0.67*)
ω=1/1.35;  (* (1 second)/(Driving period) for the video of the green marble *)
c = .01;  (* My estimate for the friction constant of the green marble *)
nT = 0; (* number of Transients: the number of driving periods integrated before plotting *)
nP =20; (* number of driving periods plotted *)
{θ0, v0} = {0, 0}; 
runPendulum2
(*nP20p1 = p1*)(* That last line was used to save the first plot *)

Here is a saved result comparing nP =20 (top) and nP = 25 (bottom).  The two runs diverge somewhere around 17 driving periods.