Shaken Washboard, or the Driven Damped Pendulum.
Version 3.1. This plots the Poincare Map.

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)];

runPoincare:= Module[{p1, p2, poincareList},
Print["Parameters are {a, c, ω} = ", {a, c, ω}];
Print["Initial conditions are {θ, θ'} = ", {θ0, v0}];
θ0 = xInIntervalab[θ0, -Pi, Pi]; (* just in case nT = 0 and θ0 is outside the interval [-π, π) *)
nT =Floor[nT];
nP = Floor[nP];  (* Just to make sure that nT and nP are integers *)
Do[
sol=NDSolve[{θ''[t]+c *θ'[t]+Sin[θ[t]]==a*ω^2*Cos[ω*t],θ[0]== θ0,θ'[0]==v0},θ,{t,0,2 Pi/ω}];
θfinal = θ[2 Pi/ω]/.sol[[1]];
θ0 = xInIntervalab[θfinal, -Pi, Pi]; (* new initial conditions after the transients *)
v0 = θ'[2 Pi/ω]/.sol[[1]] (* new initial conditions *)
, {nT}];
(*Print["The new initial conditions are the final state after nT = ", nT, " driving periods:"];*)
Print["After nT = ", nT, " transients, the final {θ, θ'} is ", {θ0, v0}];
poincareList = {};
Do[
sol=NDSolve[{θ''[t]+c *θ'[t]+Sin[θ[t]]==a*ω^2*Cos[ω*t],θ[0]== θ0,θ'[0]==v0},θ,{t,0,2Pi/ω}];
θfinal = θ[2Pi/ω]/.sol[[1]];
θ0 = xInIntervalab[θfinal, -Pi, Pi];
v0 = θ'[2Pi/ω]/.sol[[1]];
AppendTo[poincareList, {θ0, v0}]
, {nP}];
Print["After nP = ", nP," points are plotted, the final {θ, θ'} is ", {θ0, v0}];
ListPlot[poincareList, PlotRange -> All, Frame -> True, FrameLabel->{"θ", "dθ/dt"}]
]

Now, you can run the cell after defining the global variables.
The runPoincare functions does a numerical integration for every period, so nT and nP can be large without confusing Mathematica.
If you want to save the run, you can insert the cursor between the input cell and the output, and type in a comment with command 7.

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 = 300; (* number of Transients: the number of driving periods integrated before plotting *)
nP = 200; (* number of driving periods plotted (10,000 takes about 9 seconds)*)
{θ0, v0} = {-2, 1};
runPoincare

This saved run took a bit more than 10*9 seconds.  The run time scales roughly like nP

This is the youtube video parameter.
Note what happens after about 1900 driving periods.

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 = 1000; (* number of driving periods plotted *)
{θ0, v0} = {0, 0}; 
runPoincare

The conclusion is that we are observing a chaotic *transient* with the Actual Shaken Washboard.