Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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 \[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?

share|improve this question
You should get a very helpful warning about MaxSteps. – b.gatessucks Apr 26 '13 at 10:57
up vote 2 down vote accepted

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

driven damped pendulum solution

share|improve this answer
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]

enter image description here

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.