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

This site is temporarily in read only mode and not accepting new answers.

Not the answer you're looking for? Browse other questions tagged .