Solution to a T+U = E equation

Question

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.

Answer 1

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.