1
$\begingroup$

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?

$\endgroup$
1
  • 4
    $\begingroup$ You should get a very helpful warning about MaxSteps. $\endgroup$ Apr 26, 2013 at 10:57

2 Answers 2

4
$\begingroup$

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

$\endgroup$
4
  • $\begingroup$ 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$? $\endgroup$ Apr 26, 2013 at 13:22
  • 1
    $\begingroup$ @Peter, that's what ParametricNDSolve is for. $\endgroup$
    – user21
    Apr 26, 2013 at 13:29
  • $\begingroup$ Ok, looks like it is time to update to Mathematica 9 :). $\endgroup$ Apr 26, 2013 at 13:40
  • $\begingroup$ @Peter, oh yes :-) many good new NDSolve features... $\endgroup$
    – user21
    Apr 26, 2013 at 14:04
1
$\begingroup$

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]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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