This is the nature of chaos. the number of digits of precision you need to obtain convergence to time t is often exponential in t.

That said, using better integrators, it is possible to do better.

using OrdinaryDiffEq, Plots
function pend(dy, y, p, t)
    (;gamma, w, w0, B) = p
    dy[1] = y[2]
    dy[2] = gamma*w0^2*cospi(w*t) - w0^2*sin(y[1]) - 2*B*y[2]
end
# guessing at initial conditions because you didn't give them in the post
prob = ODEProblem(pend, [big"-1.0", 0], (big"0.0",50), (gamma=1.16, B=3*big(pi)/4, w0=3*big(pi), w=2))
sol1 = solve(prob, Vern7(), abstol=1e-19, reltol=1e-19);
sol2 = solve(prob, Vern7(), abstol=1e-20, reltol=1e-20);
plot(sol1, idxs=2) # blue
plot!(sol2, idxs=2) # orange

This is the nature of chaos. the number of digits of precision you need to obtain convergence to time t is often exponential in t.

That said, using better integrators, it is possible to do better.

using OrdinaryDiffEq, Plots
function pend(dy, y, p, t)
    (;gamma, w, w0, B) = p
    dy[1] = y[2]
    dy[2] = gamma*w0^2*cospi(w*t) - w0^2*sin(y[1]) - 2*B*y[2]
end
# guessing at initial conditions because you didn't give them in the post
prob = ODEProblem(pend, BigFloat[0, 0], (big"0.0",50), (gamma=1.16, B=3*big(pi)/4, w0=3*big(pi), w=2))
sol2 = solve(prob, Vern7(), abstol=1e-20, reltol=1e-20);
sol2 = solve(prob, Vern7(), abstol=1e-22, reltol=1e-22);
plot(sol1, idxs=2) # blue
plot!(sol2, idxs=2) # orange

This is the nature of chaos. the number of digits of precision you need to obtain convergence to time t is often exponential in t.

This is the nature of chaos. the number of digits of precision you need to obtain convergence to time t is often exponential in t.

That said, using better integrators, it is possible to do better.

using OrdinaryDiffEq, Plots
function pend(dy, y, p, t)
    (;gamma, w, w0, B) = p
    dy[1] = y[2]
    dy[2] = gamma*w0^2*cospi(w*t) - w0^2*sin(y[1]) - 2*B*y[2]
end
# guessing at initial conditions because you didn't give them in the post
prob = ODEProblem(pend, [big"-1.0", 0], (big"0.0",50), (gamma=1.16, B=3*big(pi)/4, w0=3*big(pi), w=2))
sol1 = solve(prob, Vern7(), abstol=1e-19, reltol=1e-19);
sol2 = solve(prob, Vern7(), abstol=1e-20, reltol=1e-20);
plot(sol1, idxs=2) # blue
plot!(sol2, idxs=2) # orange