# Derivative of function solved from NDSolve can't be plotted

I have used NDSolve to solve for a function (an angle of a triple pendulum), and now I wish to plot the derivative of that function (the angular velocity). Mathematica plots an empty plot. I had this issue before with a double pendulum system and ended up having to copy and paste cells into a new notebook one by one, and reevaluating everything. Then it suddenly started working (for the double pendulum). Doing that is not working in the case of my triple pendulum; the plot is still empty.

Needs["VariationalMethods"]
ClearAll[Global]

variables = {Subscript[θ, 1][t], Subscript[θ, 2][t], Subscript[θ, 3][t]};
Subscript[r, 1] =
Subscript[l, 1] {Sin[Subscript[θ, 1][t]], -Cos[Subscript[θ, 1][t]]};
Subscript[r, 2] = Subscript[r, 1] +
Subscript[l, 2] {Sin[Subscript[θ, 2][t]], -Cos[Subscript[θ, 2][t]]};
Subscript[r, 3] = Subscript[r, 2] +
Subscript[l, 3] {Sin[Subscript[θ, 3][t]], -Cos[Subscript[θ, 3][t]]};

Subscript[l, 1] = 1;
Subscript[l, 2] = 1;
Subscript[l, 3] = 1;
Subscript[m, 1] = 1;
Subscript[m, 2] = 1;
Subscript[m, 3] = 1;
g = 9.81; tMax = 10;
initial =
{Subscript[θ, 1] == Pi/2, Subscript[θ, 2] == Pi, Subscript[θ, 3] == Pi/2,
Subscript[θ, 1]' == Pi/2, Subscript[θ, 2]' == 0, Subscript[θ, 3]' == 0};

lagrangian = (Subscript[m, 1]/2) D[Subscript[r, 1], t].D[Subscript[r,
1], t] + (Subscript[m, 2]/2) D[Subscript[r, 2], t].D[Subscript[
r, 2], t] + (Subscript[m, 3]/2) D[Subscript[r, 3], t].D[
Subscript[r, 3], t] - g {0, 1}.(Subscript[m, 1]*Subscript[r, 1] +
Subscript[m, 2]*Subscript[r, 2] + Subscript[m, 3]*Subscript[r, 3]);
eqs = EulerEquations[lagrangian, variables, t];

sol = First[NDSolve[Join[eqs, initial], variables, {t, 0, tMax}]];
Plot[Subscript[θ, 2]'[t] /. sol, {t, 0, 10}]


and the plot is blank. Like I said, I did the exact same code for the double pendulum, except the Lagrangian was slightly different (no third mass), and I was working with fewer variables. I had the same problem; there was no plot. After copying and pasting, it started working. But no amount of copying and pasting and renaming variables is working with the triple pendulum.

Notice Subscript[θ, 2][t] /. sol /. t -> 5 gives an answer of -5.75151. But notice Subscript[θ, 2]'[t] /. sol /. t -> 5 simply spits out Derivative[Subscript[θ, 2]] (the same thing). So the derivative has not been calculated at all, it seems.

Why is there no plot?

Note: the theta function plots fine. I'm wanting to plot the derivative of each theta function. That is not working.

• Simply because Subscript[\[Theta], 2][t] doesn't match Derivative[Subscript[\[Theta], 2]][t], Change the variables in sol into Head /@ variables will fix the problem. – xzczd Nov 27 '14 at 4:11
• Wow, thank you so much. I think you should put that as an answer, because it doesn't seem to be a "small mistake" and is not intuitive for me. I still don't quite get why it doesn't work as it was before... – Sultan of Swing Nov 27 '14 at 4:15
• OK, after a second thought I retracted my closing vote, it's indeed a good example for pattern-matching. (I do suspect that the issue must have been raised before but can't find one.) It's a little painful for me to explain this issue systematically, let's wait for a moment first to see if any pattern-matching master will come to give an answer. – xzczd Nov 27 '14 at 6:11
• @xzczd I mention the issue in my answer here but do not give an extensive explanation. I'm pretty sure it has come up a couple of times before, too. – Michael E2 Nov 27 '14 at 14:03

Plot[D[variables /. sol, t] /. t -> u, {u, 0, 10}, Evaluated -> True] • Maybe you can go a little deeper on this issue so we will be able to mark those subsequent similar question as a duplicate of this confidently? – xzczd Nov 28 '14 at 14:48

Change the variables in sol into Head /@ variables.

Before answering your question, I'd like to suggest you to try to simplify your code sample as much as possible before posting it here so your question will draw more attention. To summarize, something seems to be "wrong" when you try to plot the derivative of the output of the NDSolve, then can we reproduce the error in a simpler case? Of course:

sol = NDSolve[y'[x] == x && y == 1, y[x], {x, 0, 1}]
(* {{y[x] -> InterpolatingFunction[{{0., 1.}}, <>][x]}} *)
Plot[y'[x] /. sol, {x, 0, 1}]
(* Blank plot *)


The sol is nothing but a rule with the form y[x] -> …[x], then is InterpolatingFunction the root of the problem? Of course not:

sol = y[x] -> Sin[x]
Plot[y'[x] /. sol, {x, 0, 1}]
(* Blank plot *)


So the question boils down to why ReplaceAll(/.) doesn't use y[x] -> … to rewrite y'[x].

To answer the question, see the following examples first:

Clear[a, y, x]
1 /. 1 -> a
(* a *)
1 + 1 /. 1 -> a
(* 2 *)


The first replacement successes, while the second "fails", why? Because the replacement done by ReplaceAll (and all the other pattern-matching functions) is literal rather than semantic i.e. things to be replaced must exist in the right hand side exactly, to be precise, exist in the FullForm:

FullForm[y'[x]]
(* Derivative[y][x] *)


Apparently y[x] doesn't exist in y'[x].

So how to fix the problem?

One way is to do the replacement before derivation, as belisarius does in his answer.

Another way is to make use of the other syntax of the 2nd argument of NDSolve: y[x] isn't necessary, a single y is enough:

sol = NDSolve[y'[x] == x && y == 1, y, {x, 0, 1}]
(* {{y -> InterpolatingFunction[{{0., 1.}}, <>]}} *)
Plot[y'[x] /. sol, {x, 0, 1}]
(* Desired plot *)


This is what I've done in the comment above: Head /@ {y[x]} gives {y}.

Actually, there exists a 3rd solution: the 2nd argument of NDSolve owns a not-well-known form, it can output the derivation directly:

sol = NDSolve[y'[x] == x && y == 1, y'[x], {x, 0, 1}]
(*{{y'[x] -> InterpolatingFunction[{{0.,1.}}, <>][x]}}*)
Plot[y'[x] /. sol, {x, 0, 1}]
(*Desired plot*)


If you still feel confused, please copy all the code into Mathematica and stroke F1 to check the document of the related document.

Though named as Long Answer, the above introduction still only covers a small aspect of pattern-matching, if you want to go even deeper, searching in this site or have a look at Leonid's excellent book.

Works for me if I change your plot line to

Plot[variables /. sol, {t, 0, 10}]
`
• I want to plot the derivative of the theta function(s), not the theta functions themselves. – Sultan of Swing Nov 27 '14 at 2:56
• Then: derivs = D[variables /. sol, t]; Plot[derivs, {t, 0, 10}] – Bill Nov 27 '14 at 6:08