# 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 β 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?

• You should get a very helpful warning about MaxSteps. Apr 26, 2013 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$? Apr 26, 2013 at 13:22
• @Peter, that's what ParametricNDSolve is for.
– user21
Apr 26, 2013 at 13:29
• Ok, looks like it is time to update to Mathematica 9 :). Apr 26, 2013 at 13:40
• @Peter, oh yes :-) many good new NDSolve features...
– user21
Apr 26, 2013 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] == 0,
x'[0] == 0}, x, {t, 0, tend}, MaxSteps -> Infinity]


Plotting a smaller subrange helps:

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