# Evaluating the solution of a differential equation

n = 3;
g = 9.8;
m = {};
l = {};
EOM = Array[# &, n];
For[i = 1, i <= n, ++i,
AppendTo[m, 0.1];
AppendTo[l, 1];]

\[Theta][t_] := Array[Subscript[\[Theta], #][t] &, n];

L[t_] = 0.5*Sum[(n + 1 - i)* m[[i]]*(l[[i]]^2* (Subscript[\[Theta], i]'[t])^2 + Sum[2*l[[i]]*l[[k]]*Subscript[\[Theta], i]'[t] * Subscript[\[Theta], k]'[t] * Cos[Subscript[\[Theta], k][t] - Subscript[\[Theta], i][t]], {k, 1, i - 1}]), {i, 1, n}] + g*Sum[  (n + 1 - i)*m[[i]]*l[[i]]*Cos[Subscript[\[Theta], i][t]]   , {i, 1, n}];

For[i = 1, i <= n, ++i,
EOM[[i]] = D[D[L[t], Subscript[\[Theta], i]'[t]], t] - D[L[t], Subscript[\[Theta], i][t]] == 0]

thetainicial = Array[# &, n];
Dthetainicial = Array[# &, n];
For[i = 1, i <= n, ++i,
thetainicial[[i]] = Subscript[\[Theta], i][0] == Pi/4;
Dthetainicial[[i]] = Subscript[\[Theta], i]'[0] == 0]
EQ = Join[EOM, thetainicial, Dthetainicial];

s = Flatten@NDSolve[EQ, \[Theta][t], {t, 0, 100}, Method -> {"EquationSimplification" -> "Residual"} ]

p[t_] := Array[Subscript[p, #][t] &, n]
For[i = 1, i <= n, ++i,
Subscript[p, i][t] = D[L[t], D[Subscript[\[Theta], i][t], t]]]


I need to do a Plot p[t] and a ParametricPlot for p[t] and [Theta][t], I'm trying like this but it isn't working:

Plot[Evaluate[p[t]/.s], {t, 0, 10}, PlotRange -> All, PlotLegends -> Automatic]

ParametricPlot[ Evaluate[p[t]/.s,\[Theta][t] /. s ], {t, 0, 10}, AspectRatio -> 1, PlotLegends -> Automatic]


The last one makes this Error: ParametricPlot::plln: Limiting value InterpolatingFunction[{{0.,100.}},<<3>>,{Automatic}][t] in {<<1>>} is not a machine-sized real number.

Can anyone help me?

• It would be best to rewrite the code without using subscripted variables. They create problems. – David Keith Jun 24 '17 at 18:38
• Not really sure how to do that – ines Jun 24 '17 at 19:06
• I often use for example n[1] as a subscripted n. These work more reliably as variables. – David Keith Jun 24 '17 at 20:05

Another approach to obtaining the derivatives of θ is

s = Flatten@NDSolve[EQ, Flatten@{θ[t], D[θ[t], t]}, {t, 0, 100},
Method -> {"EquationSimplification" -> "Residual"} ]


The basic problem is that s does not contain a replacement rule for the derivatives. Then for the parametric plot, there is also the problem that the first argument is not of the form {{fx,fy}, {gx,gy}...}. Fixing this, we get:

rules = Join[s, Table[Derivative[1][Subscript[\[Theta], i]][t] ->
D[s[[i, 2]], t], {i, 3}]];
Plot[Evaluate[p[t] /. rules], {t, 0, 10}, PlotRange -> All,
PlotLegends -> Automatic]


and

fns = Transpose[{{Subscript[\[Theta], 1][t],
Subscript[\[Theta], 2][t], Subscript[\[Theta], 3][t]}, p[t] /. rules}];
ParametricPlot[Evaluate[fns /. rules], {t, 0, 10}, AspectRatio -> 1,
PlotLegends -> Automatic]


The following line deserves a bit of extra attention:

Table[Derivative[1][Subscript[\[Theta], i]][t] -> D[s[[i, 2]], t], {i, 3}]]


Why did I write the replacement rule like this? The way I did it was that I evaluated p[t]:

and then copied one of the derivatives and asked for the input form:

then I knew what the left-hand side of the replacement rule should be.

• @ines In principle, this should be easy to scale. In fns, Table will have to be used to create the required number of functions. Where I have written 3 in some places, that will have to be replaced with n. Also make sure to read bbgodfrey's answer and understand it, it can help make things easier for you. – C. E. Jun 24 '17 at 22:36
• Thank you for your kind comment. – bbgodfrey Jun 25 '17 at 4:02