How to improve the solution of a double pendulum system

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}]


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.

-
The problem is with plotting, not NDSolve. Try setting the options PlotPoints -> 100, MaxRecursion -> 8 on ParametricPlot. – Michael E2 Nov 18 '13 at 14:04
Right, this solves it. Any idea why it doesn't plot smooth functions by default? – shrx Nov 18 '13 at 14:06
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. – Michael E2 Nov 18 '13 at 14:34
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. – Szabolcs Nov 18 '13 at 18:26

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]


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]


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.

-
@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. – Michael E2 Nov 10 '14 at 13:36
@Phonon If you like, but, quite sincerely, your gratitude is enough of a reward. :) – Michael E2 Nov 10 '14 at 13:42