Hello, how to build a phase diagram (p versus q) with the graphs that are generated in the commands? Identical, for example, to the ones that appear when we use the command: ParametricPlot[{Sin[u], Sin[2 u]}, {u, 0, 2 Pi}]
H = p^2/2 + (1 - Cos[q]);
a = 0;
b = 20;
p[0] = 1;
h = 0.5;
L = (b - a)/h;
Do[q[n] = n*h, {n, 0, L}]
f[q_, p_] = D[H, p];
Do[{
k1 = h*f[q[n], p[n]],
k2 = h*f[q[n] + h/2, p[n] + k1/2],
k3 = h*f[q[n] + h/2, p[n] + k2/2],
k4 = h*f[q[n] + h, p[n] + k3],
p[n + 1] = p[n] + (k1/6 + k2/3 + k3/3 + k4/6)},
{n, 0, L}];
DiscretePlot[{p[n]}, {n, 0, 20}, Joined -> True, Filling -> None,
AxesLabel -> {q, p}]
Do[q[n] = n*h, {n, 0, L}]
g[q_, p_] = -D[H, q];
Do[{
k1 = h*f[q[n], p[n]],
k2 = h*f[q[n] + h/2, p[n] + k1/2],
k3 = h*f[q[n] + h/2, p[n] + k2/2],
k4 = h*f[q[n] + h, p[n] + k3],
p[n + 1] = p[n] + (k1/6 + k2/3 + k3/3 + k4/6)},
{n, 0, L}];
DiscretePlot[{q[n]}, {n, 0, 20}, Joined -> True, Filling -> None,
AxesLabel -> {q, p}]