NDSolve solution for driven damped pendulum diverges

Question

I want to solve numerically for the system of the driven damped pendulum using Mathematica. This is the second-order nonlinear equation \begin{equation} \ddot{x} + 2 \beta \dot{x}+ \omega_0^2 \sin x = \gamma \omega_0^2 \cos[\omega t] \end{equation} The Mathematica code I used is simply

s = NDSolve[{x''[t] + 2 β x'[t] + ω0^2  Sin[x[t]] == γ  ω0^2 Cos[ω t], x[0] == 0, x'[0] == 0}, x, {t, 0, 500}, PrecisionGoal -> 20][[1]];

With parameters

β = ω0/4; γ = 0.2; ω = 2 π; ω0 =  1.5 ω; 

However, the solution

x[t] /. s 

Is an interpolation function defined only for $t\in [0,101]$ roughly. Plotting the function shows that it suddenly diverges at $t \approx 101$, which clearly should not happen. I tried changing the PrecisionGoal but this only change the exact value at which the function diverges - it still diverges at some value well before $t=500$. The problem persists for other values of the parameters, when I change $\gamma$.

Is there any way to get an accurate solution for longer times?

Answer 1

Using exact values (e.g. setting γ = 1/5 instead of γ = 0.2) and increasing the MaxSteps option should yield a reasonable result. Below, I add the filigree of also using a method I tend to prefer for problems with oscillatory solutions (although I've found that just setting Method -> "StiffnessSwitching" without increasing MaxSteps still works well):

xp = With[{γ = 1/5, ω = 2 π},
     Block[{ω0 = 3 ω/2, β}, β = ω0/4;
           x /. First @ 
           NDSolve[{x''[t] + 2 β x'[t] + ω0^2 Sin[x[t]] == γ ω0^2 Cos[ω t], 
                    x[0] == x'[0] == 0}, x, {t, 0, 500}, MaxSteps -> 1*^5, 
                   Method -> "StiffnessSwitching", PrecisionGoal -> 20]]]

Inspect a few sections:

{{Plot[xp[t], {t, 0, 50}], Plot[xp[t], {t, 100, 150}]},
 {Plot[xp[t], {t, 350, 400}], Plot[xp[t], {t, 450, 500}]}} // GraphicsGrid

Answer 2

The issue is not with NDSolve but plotting:

β=ω0/4;γ=1/5;ω=2 π;ω0=3/2 ω;
 tend = 500;
s = NDSolveValue[{x''[t] + 
     2 β x'[t] + ω0^2 Sin[
       x[t]] == γ ω0^2 Cos[ω t], x[0] == 0, 
   x'[0] == 0}, x, {t, 0, tend}, MaxSteps -> Infinity]

Plotting a smaller subrange helps:

Plot[s[t], {t, 300, 400}, PlotRange -> All]