# Solution to a T+U = E equation

I needed to solve really easy differential equation (in dimensionless units): $$\mathcal{T} (\dot{\xi}) + \mathcal{U} (\xi) = \text{const.} \equiv \varepsilon ; \quad T(\dot{\xi}) = \dot{\xi}^2 ; \quad \left( \quad \right)^{\dot{\;}} \equiv \frac{\mathrm{d}}{\mathrm{d} \tau}$$ and $\mathcal{U} (\xi)$ is guaranteed to have a local minimum somewhere (possibly more than one).

I wanted to keep it as simple as I could, so the code looks something I would write down at the paper, the only help I needed from Matematica is the one with non-elementary functions:

U[ξ_] := ξ^2; (*1*)

f[ξ_, Subscript[x, 0] _, E_] :=
NIntegrate[1/Sqrt[E - U[u]], {u, Subscript[x, 0], ξ}]; (*2*)

g[Subscript[ξ, 0], E_, t_] = InverseFunction[f, 1, 3]; (*3*)

Plot[g[1, 2, t], {t, 0, 2}] (*4*)


I've numbered the lines for easier reference. The first line (1) defines our potential energy to be a quadratic potential (however, some more complex one could be chosen, in fact, quadratic potential is the most simple one and solvable in elementary functions). The law of conservation of energy states: $$\dot{\xi}^2 + \mathcal{U} (\xi) = \varepsilon \implies \int_{\xi_0}^{\xi} \frac{\mathrm{d} u}{\sqrt{\varepsilon - \mathcal{U} (u)}} = \tau$$

So, basically, second line (2) defines function $f$: $$f(\xi, \xi_0, \varepsilon) \equiv \int_{\xi_0}^{\xi} \frac{\mathrm{d} u}{\sqrt{\varepsilon - \mathcal{U} (u)}}$$

I thought that for general potential it's impossible to find primitive function, so it's not ever worth a try to write Integrate, so I used NIntegrate instead.

After this step we might have the solution in implicit form, but the explicit form would be more appropriate: $$f(\xi, \xi_0, \varepsilon) = \tau \implies \xi = g( \xi_0, \varepsilon, \tau)$$

That's what I've tried to do in third line (3). The documentation states, that "InverseFunction[f,n,tot] represents the inverse with respect to the $n^{(th)}$ argument when there are $tot$ arguments in all.". Final step is to plot this solution at line (4) for some initial $\xi_0 = 1$, $E = 2$ aaaaand...what a surprise, empty plot!

Now I really love Mathematica, but things like this drives me crazy. What should I do? You know, if I wanted to solve motion in quadratic potential I wouldn't need Mathematica, this is just a test potential, the goal is to obtain solution for some more complex potential (biquadratic, cosine, ...). I have the suspicion, that the whole problem can be solved automatically by the fancy instructions like Solve, DSolve, NDSolve etc., but I don't understand those complex instructions like Plot[f/.solution] or other strange symbols in code (&, \$, //, /., /@, #, ...), neither can I find them in the documentation. Especially /. is everywhere (like something /. something else - and this magically does the right thing), but no one ever explains what it does.

Any help would be greatly appreciated.

• /. is a shorthand for ReplaceAll[]. You may try to read this question & answers: mathematica.stackexchange.com/q/18393/193 – Dr. belisarius Sep 18 '14 at 1:29
• Unfortunately, this requires a long, long answer which I have no energy to give. But, I'd first like to advise a deep breath as everyone encounters this stage when learning Mathematica. It does make sense at some point. Next, read though the answers to this question. After that, go through the doc examples, execute them as is, and modify them until they break. Then, fix the breaks. – rcollyer Sep 18 '14 at 1:31
• For /., this answer in particular: Using the result of functions that return replacement rules. – Rahul Sep 18 '14 at 1:34
• @rcollyer +14 for recommending a deep breath. We all have been there. The other alternative is crying out loud. Not recommended at all. – Dr. belisarius Sep 18 '14 at 2:03
• To all: deep breath and crying out loud, that's what I'm being through these days considering how difficult is Mathematica's syntax to understand for me. At least you made me grin...at least I know these are the problems everyone encounter at the beginning. – user16320 Sep 18 '14 at 8:52

Your method is perfectly fine in principle, because in one dimension energy conservation leads to a "quadrature" solution (i.e., you formally avoid solving the usual equation of motion). But since inverting the numerical integration is another numerical step it doesn't actually end up making the solution very easy to obtain. For example, you have to think about the range of initial conditions consistent with your chosen energy first.

In more than one dimension, this approach based solely on energy conservation is not sufficient to pin down the solution. So I would suggest a more generalizable approach that starts from the two terms (kinetic and potential energy) and automatically derives the correct equations of motions without directly making use of energy conservation as you do in your question.

One approach that does this is the Lagrangian method. Here is how it can be easily programmed in Mathematica:

Needs["VariationalMethods"]

ClearAll[ξ, t];
U = ξ[t]^2;
T = D[ξ[t], t]^2;

solution =
First@NDSolve[
Flatten[
{EulerEquations[T - U, ξ[t], t], ξ[0] == 0, ξ'[0] == 1}
], ξ[t], {t, 0, 10}];

Plot[ξ[t] /. solution, {t, 0, 10}]


Here I give NDSolve the equations of motion as calculated by EulerEquations, together with a set of initial conditions for the position ξ[0] == 0, and velocity ξ'[0] == 1. What goes into EulerEquations is the Lagrangian T - U, the function ξ[t] describing the trajectory you want to solve for, and the name of the time variable t. The Flatten command is only needed when you generalize to higher dimensions where there will be more equations of motion and more initial conditions. Since NDSolve and similar function in general output a list of lists {{...},...{...}}, I added the command First@NDSolve to select the first sublist which contains a single complete solution. Having the result wrapped in an extra list can sometimes lead to unexpected results in later manipulations.

It's also quite straightforward to implement Hamilton's equations of motion if you want to write the trajectory in terms of position and momenta in phase space. But your question seems to indicate that you are not doing that.

Using InverseFunction

If you want to compare the above to your original approach, I would do it as follows:

U[ξ_] := ξ^2;

Clear[f];
f[ξ_?NumericQ, ξ0_, e_] :=
NIntegrate[1/Sqrt[e - U[u]], {u, ξ0, ξ}];

Clear[g];
g[pos_, speed_] := InverseFunction[f[#, pos, U[pos] + speed^2] &]

trajectory[.1]

(* ==> 0.0998334 *)

trajectory[1.5]

(* ==> 0.997495 *)


Here, I defined the first integral to be evaluated only when the argument ξ is numerical, so that InverseFunction knows not to try and symbolically extract the inverse (which would fail because the whole function is based on a numerical integration).

In the function g that uses InverseFunction, the two arguments are the initial position and the absolute value of the velocity. They are then used to calculate the appropriate value of the energy e. Notice that I replaced the capital E by e. This is one of the first things you should do:

• Don't use capitalized variable names unless you know exactly what you're doing (or you work in marketing and hence don't care what you're doing).
• Don't use Subscript for variable names (same exceptions as above)

The capital E is a reserved symbol. If you do want capitals, enter the E as EscapeEEscape instead.

Now g is itself a function of a single argument, which I test by invoking it with an extra [t] at the end. To make the distinction clearer, I give the result of g[0,1] a name, trajectory, and then call trajectory[t] to get an actual value.

However, this inverse function approach breaks down as soon as you hit the first classical turning point. So if you try to use this method for times larger than the first traversal time to the turning point, it will not work. You see this in the last calculation in the significant slowdown for t=1.5. This is why I suggested an alternate approach earlier.

• I'm very grateful for this short code that works, but you've just fed the hungry man. I'd be even more grateful if you teach the man how to hunt, so: why do you need to add command First@ before NDSolve? I've tried to play a little with more potentials (like pendulum U = cos(x) potential) and I deleted First@ and it works anyway. In the documentation I found that First[expr] gives the first element in expr. How is this useful in here? Thanks, your answer is great, although I don't fully understand it. – user16320 Sep 18 '14 at 8:58
• I've added some more explanation and also addressed your original InverseFunction` approach in my edited answer. – Jens Sep 18 '14 at 18:46
• I'm really grateful for you patience, this answer is so much helpful...Can't thank you enough! Now I've managed to produce working code for some animation of oscillators in phase space diagrams in time for various interesting potentials and even add dissipation term to demonstrate when the Liouville theorem will work and when it won't. – user16320 Sep 18 '14 at 23:10