# 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 \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, x' == 0}, x, {t, 0, 500}, PrecisionGoal -> 20][];


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?

• You should get a very helpful warning about MaxSteps. – b.gates.you.know.what 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 == x' == 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:

β=ω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,
x' == 0}, x, {t, 0, tend}, MaxSteps -> Infinity]


Plotting a smaller subrange helps:

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