# NDSolve solution for driven damped pendulum diverges

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

s = NDSolve[{x''[t] + 2 \[Beta] x'[t] + \[Omega]0^2  Sin[x[t]] == \[Gamma]  \[Omega]0^2 Cos[\[Omega] t], x[0] == 0, x'[0] == 0}, x, {t, 0, 500}, PrecisionGoal -> 20][[1]];


With parameters

\[Beta] = \[Omega]0/4; \[Gamma] = 0.2; \[Omega] = 2 \[Pi]; \[Omega]0 =  1.5 \[Omega];


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?

-
You should get a very helpful warning about MaxSteps. – b.gatessucks Apr 26 '13 at 10:57

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


-
Thanks! I now want to create a list of solutions for different values of $\gamma$. Simply defining a table with the command above and $g$ replacing 1/5, running from, say, 1.06 to 1.08 in steps of 0.001 doesn't work. Is there any way I can create a list of solutions for different $\gamma$? – Yiteng Apr 26 '13 at 13:22
@Peter, that's what ParametricNDSolve is for. – user21 Apr 26 '13 at 13:29
Ok, looks like it is time to update to Mathematica 9 :). – Yiteng Apr 26 '13 at 13:40
@Peter, oh yes :-) many good new NDSolve features... – user21 Apr 26 '13 at 14:04

The issue is not with NDSolve but plotting:

\[Beta]=\[Omega]0/4;\[Gamma]=1/5;\[Omega]=2 \[Pi];\[Omega]0=3/2 \[Omega];
tend = 500;
s = NDSolveValue[{x''[t] +
2 \[Beta] x'[t] + \[Omega]0^2 Sin[
x[t]] == \[Gamma] \[Omega]0^2 Cos[\[Omega] 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]


-