4
$\begingroup$

First, the code:

(* numerical solution of the double pendulum system in polar coordinates *) sol = Flatten[{v1, v2, θ1, θ2} /. NDSolve[{v1'[t] == -l^2/ 2 (3 g/l Sin[θ1[t]] + θ1'[t] θ2'[ t] Sin[θ1[t] - θ2[t]]), v2'[t] == -l^2/ 2 (g/l Sin[θ2[t]] + θ1'[t] θ2'[ t] Sin[θ1[t] - θ2[t]]), θ1'[t] == 6/(m l^2) ( 2 m v1[t] - 3 Cos[θ1[t] - θ2[t]] m v2[t])/( 16 - 9 Cos[θ1[t] - θ2[t]]^2), θ2'[t] == 6/(m l^2) ( 8 m v2[t] - 3 Cos[θ1[t] - θ2[t]] m v1[t])/( 16 - 9 Cos[θ1[t] - θ2[t]]^2), θ1[0] == π/2, θ2[0] == π/2, v1[0] == v2[0] == 0 } /. {g -> 9.81, l -> 1, m -> 1}, {v1, v2, θ1, θ2}, {t, 0, 20} ]] (* cartesian coordinates of pendulum ends *) x1k[t_, l_] := l Sin[sol[[3]][t]] y1k[t_, l_] := -l Cos[sol[[3]][t]] x2k[t_, l_] := l (Sin[sol[[3]][t]] + Sin[sol[[4]][t]]) y2k[t_, l_] := -l (Cos[sol[[3]][t]] + Cos[sol[[4]][t]]) (* corresponding plot *) ParametricPlot[ Evaluate[{{x1k[t, l], y1k[t, l]}, {x2k[t, l], y2k[t, l]}} /. {l -> 1}], {t, 0, 20}]

double pendulum plot

As you can see, the function for the position of the second pendulum's end is very inaccurate, so the red curve is jagged instead of smooth. I have tried with specifying various options to NDSolve: AccuracyGoal, InterpolationOrder, MaxStepFraction, PrecisionGoal, but neither seemed to help.

Improve this question
$\endgroup$
4
  • 3
    $\begingroup$ The problem is with plotting, not NDSolve. Try setting the options PlotPoints -> 100, MaxRecursion -> 8 on ParametricPlot. $\endgroup$ Commented Nov 18, 2013 at 14:04
  • $\begingroup$ Right, this solves it. Any idea why it doesn't plot smooth functions by default? $\endgroup$ Commented Nov 18, 2013 at 14:06
  • $\begingroup$ Because MaxRecursion limits how much effort the plotting function may put into making it smooth. Some functions would take a very long time to make smooth. $\endgroup$ Commented Nov 18, 2013 at 14:34
  • 1
    $\begingroup$ The number of points that will be plotted is limited for performance reasons, as in all other systems. Mathematica actually has an advantage here because it does adaptive plotting. First it divides the {t,0,20} interval into PlotPoints number of equal parts and computes the function at each point. Many other systems will do only this much, and will typically give you some jaggedness. But Mathematica will look at each segment pair and if their angle is large enough (i.e. if they look jagged), it will subdivide them to make the plot smoother. Then it repeats this up to MaxRecursion times. $\endgroup$ Commented Nov 18, 2013 at 18:26

1 Answer 1

Reset to default
6
+50
$\begingroup$

The problem has to do with the sample points used to contruct the plot, and not with NDSolve. Mathematica automatically subdivides segments when the angle is greater than some limit, but it will do so only MaxRecursion times. One can increase PlotPoints to increase the number of initial sample points and increase MaxRecursion to let Mathematica subdivide further where needed. 

ParametricPlot[ Evaluate[{{x1k[t, l], y1k[t, l]}, {x2k[t, l], y2k[t, l]}} /. {l -> 1}], {t, 0, 20}, PlotPoints -> 100, MaxRecursion -> 8]

Mathematica graphics

The angles one sees in the plot in the question usually result from a large velocity, or in other words from a large change in arclength for a given change in parameter t.

One can see that is the case here from the plot of the speed:

Plot[Evaluate @ Norm @ {D[x2k[t, 1], t], D[y2k[t, 1], t]}, {t, 0, 20}, PlotRange -> All]

Mathematica graphics

Another approach, perhaps not needed here, is to reparametrize by arclength. Then the velocity will be constant. If the curve is not very long, then the default plotting options will probably produce a smooth graph.

For this problem the speed is zero at t == 0, which produces a singularity if we use the standard initial value t[0] == 0. So instead we integrate the speed over a short interval to get a nonzero initial value.

invarclength = Quiet @ NDSolveValue[ {t'[s] == 1/Norm @ {Derivative[1, 0][x2k][t[s], 1], Derivative[1, 0][y2k][t[s], 1]}, t[NIntegrate[Norm @ {Derivative[1, 0][x2k][t, 1], Derivative[1, 0][y2k][t, 1]}, {t, 0, 10^-16}]] == 10^-16, WhenEvent[t[s] == 20, "StopIntegration"]}, t, {s, 0, Infinity}] ParametricPlot[ Evaluate[{{x1k[#, 1], y1k[#, 1]}, {x2k[#, 1], y2k[#, 1]}} &@ invarclength[s]], Evaluate @ {s, Sequence @@ First @ invarclength["Domain"]}]

It produces a plot similar to the one above, but of about 1/3 the ByteCount.

Improve this answer
$\endgroup$
2
  • $\begingroup$ @Phonon If ImageSize is Automatic, then image sizes are automatically scaled in certain contexts. You could use GraphicsRow[{##}] & instead of Row. Or use the ImageSize option. Or both. $\endgroup$ Commented Nov 10, 2014 at 13:36
  • 2
    $\begingroup$ @Phonon If you like, but, quite sincerely, your gratitude is enough of a reward. :) $\endgroup$ Commented Nov 10, 2014 at 13:42

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.