Solving equation of motion for position dependent force

The plan here is to set up step by step the problem in terms of Hamiltonian mechanics. The principal idea are that the Hamiltonian $H$ represents the energy of the system $$H = T + U \quad \text{(kinetic + potential energy)}$$ and the differential equations governing the motion are derived from $H$ as follows: $$ dp/dt = - \partial H/\partial x, \quad dx/dt = \partial H/\partial p$$ where $x$ and $p$ are the position and momentum variables respectively. In the problem at hand the kinetic energy is given by $T = p^2/(2m)$.

One feature of Hamiltonian systems is that $H$ is invariant. It is a so-called "first integral" (which just means is constant as the system evolves). And when the system comes from a physical system, its invariance is equivalent to the conservation of energy. Consequently, the solution to an initial value problem $x(t_0)=x_0,\ p(t_0)=p_0$ traces the curve $$H(x,p) = H_0 \buildrel {\rm DEF} \over = H(x_0,p_0)$$ in phase space. If the curve is closed, then the solution must be periodic (assuming $H$ is nice enough).

U[x_] := x^-2 - 2/x; (* normally single/initial capital letters are *) H[x_, p_] := p^2/(2 m) + U[x]; (* to be avoided, but I'm going to break that rule *) Plot[U[x], {x, 0, 7}, GridLines -> {None, {H[0.55, 0.]}}] 

An initial condition that gives a negative value for $H_0$ leads to periodic solution because the total energy $H=p^2/(2m)+U$ cannot be less than the potential energy $U$:

H[0.55, 0.] (* -0.330579 *) Block[{m = 0.01}, ContourPlot[H[x, p] == H[0.55, 0.], {x, 0, 6}, {p, -0.2, 0.2}, FrameLabel -> Automatic, PlotLabel -> HoldForm[H == Style[H0, Italic]]]] 

One can calculate the period by integrating $dt$ over the curve, where in the Hamiltonian system we have $$dt = {dx \over dx/dt} = {1 \over \partial H/\partial p} \; dx = {m \over p}\;dx$$ In our case, if $\gamma$ is the curve $H = H_0$, we have $$\int_\gamma {1 \over \partial H / \partial p} \; dx = 2 \int_{x_{\rm min}}^{x_{\rm max}} {m \over \sqrt{2 m \, (H_0 - U(x))}} \; dx$$ If we solve the differential equation over one period, then we extend the solution to all time by periodicity.

All this can be done in Mathematica.

• Input: u, m, x0, v0 or p0. I Rationalize[] the numeric parameters because I used some exact solvers (Solve[], Integrate[]). Numeric solvers NSolve[], NIntegrate[] could be substituted.
• Construct Hamiltonian ham.
• Set up the IVP system (ODE & ICs). I found it convenient to solve these equations for their principal variables. This let's us substitute the "right-hand sides" of the system into expressions like the Hamiltonian and period integral.
• Calculate the invariant value of the Hamiltonian ham0 ($h_0$).
• Solve ham0 - u == 0 for the range of $x$, {xmin, xmax}`. For a general solver, this step would need some fixing. It's possible there is more than one loop, that the curve is unbounded, or consists of a single point.
• Construct dt in terms of x.
• Integrate dt to get period.
• Solve the IVP with NDSolve[] over {t, 0, period}.
• The periodic solution can be code x[Mod[t, period]] /. sol. One potential issue is that the period will be perturbed by the numeric method in NDSolve[]. It's not a problem in this case.

I've commented out some localizations so one can examine what the results are after executing. One can localize some or all them as desired. One could make a function that would take u, m, x0, v0 as inputs, provided one decides some user-interface questions.

Module[{m, x0, v0, u, f(*,rhsSub,xmin,xmax,sys,ham,h0,period,dt,x,p*)}, (* Inputs *) u = U[x[t]]; (* potential (given) *) {m, x0, v0} = Rationalize@{0.01, 0.55, 0}; (* parameters/conditions (given) *) (* Set-up *) ham = u + p[t]^2/(2 m); (* Hamiltonian *) f = -D[u, x[t]]; (* force -- unnecessary *) sys = { (* ivp system *) p'[t] == -D[ham, x[t]], (* force does appear here, too *) x'[t] == D[ham, p[t]], x[0] == x0, p[0] == m*v0}; (* initial conditions *) rhsSub = First@Solve[sys, {x'[t], p'[t], x[0], p[0]}]; (* solution rules for sys *) ham0 = ham /. t -> 0 /. rhsSub; (* initial value of the Hamiltonian *) {xmin, xmax} = x /. Solve[ (* range of x on (closed) trajectory in phase space *) ham0 - u == 0 /. x[t] -> x]; dt = (* differential of time *) 1/x'[t](* times dx*)/. rhsSub /. {p[t] -> Sqrt[2 m (ham0 - u)]}; period = 2 Integrate[dt /. x[t] -> x, {x, xmin, xmax}]; {ndsol} = NDSolve[sys, {x, p}, {t, 0, period} (*opts*), InterpolationOrder -> All] (* Use InterpolationOrder -> All if you want more accuracy in the InterpolatingFunction[] *between* the steps *) ] 

A plot over three periods:

Plot[x[Mod[t, period]] /. ndsol, {t, 0, 3 period}] 

Another view of the evolution of the system: the position X (x[Mod[t, period]]), a plot of the potential energy U at time t (U[x[Mod[t, period]]]), and a plot of the point {x[Mod[t, period]], p[Mod[t, period]]} tracing out the integral curve H == H0 (ham == ham0).

Code image (Manipulate[] version)

Uncompress@"1:eJzlV81uEzEQTlrSQqnopRckDvAGvEEIgbaREoiyqeCIs+vNWvWuF9urJu2Vl4A34MIrcEaCC2/QR+CKkMBjZ+PsZjdBIfxE+DCxxzPj8cw3s869Aev51UqlIvYVOWHUa7Iwplhifwu41xR5RmRgV20ipFltK+JgKWrqN4gpm7Bh+ZAy98zYLdaBrZBfHr47vrp4Vfd3Ye+mIk0WSZbwbsba45cJomYJeieajjSNZ9j8zWsYV3VemQx/e9YDo0Zgg+wW78YTGx/ep5O6vXkvoVjcUJNWiIbYIReYfAUzWYH9VKCLPI9Ew6LIpStyG3zZyhjR/APg1+aNw+kQmx6KhniRZXuPTx9h3HlQtl9yT0hGQ8TYVWdJwmwy+iTEglSnvgGvy84xJ6BLvqsxD49Ew0Nn+RbII45CLDlxbaKzWTLYAXI6YRSkkn97++VzZ9Ctm3TOh+qYE69NIixyR8DkCYtwEULLoLTu8PxhRNXWjyit20gkC1UIXCO9p0gHRSROKEp7yI6ua5qEUTaJPgiL6zpPKA6IK3JJhlUPeyZL2lFGIglh0mFtIz6cpBDy/2gcoVD5UbVhV9LWrR6OKXJxg9LcNadQG1ngdpinWVJ7H6v0MU93Dr0rGM2FESaNERZ60ucJ/r1hntqFex5BMZUdfJAK9Il7Jhw5pguPt+F8GiOXyLG+dAk3W1Sm7yN1ynSW82Uv48tyNwDLzQBxqUD/wnERxV6ZuoZRy8ORVN5lFus0tVE+/gPNZYUuCRgp65JQQG00wNTEC64NL5YjxkNdrM/9nVTHCdj5zJdnKr8ZnWbm6/fTOiWNaZUGsbYCPpxH9N+vi42q3WWQP/Vrc5APNhvy8Uo4/49h/KvuLXxgLAFh7k+aeT1bpskDLE/u6xy2JKLEdUDsbsGTXs48LOx3QstZBPwAu4ZrPw==" 

Here is a check on how the Hamiltonian (energy) is conserved. (This is the reason for InterpolationOrder -> All. Without it, there would be greater interpolation error between the steps.)

Plot[ham - ham0 /. t -> Mod[t, period] /. ndsol // RealExponent // Evaluate, {t, 0, 2 period}, PlotPoints -> 201, PlotRange -> {-16, 0}, GridLines -> {None, {-7, -8}}] 

Other available tools

1. When a system has an invariant, one can use the projection method in NDSolve[]. In our case the option code would be something like either of these:

Method -> {"Projection", "Invariants" -> {ham}} Method -> {"Projection", "Invariants" -> {ham}, Method -> "ExplicitRungeKutta"} 

After each step, the solution is adjusted by projecting it onto curve ham == <constant>, where the constant is calculated by NDSolve[] during initialization (it will be equal to ham0). The first one actually has a couple of glitches; one can try another method, as shown in the second option. The plot of the invariant for "ExplicitRungeKutta" is shown below. It is quite a bit better.

Plot[ham - ham0 /. t -> Mod[t, period] /. ndsol // RealExponent // Evaluate, {t, 0, period}, PlotPoints -> 201, PlotRange -> {-18, 0}, GridLines -> {None, {-7, -8}}] 

2. Another method designed specially for Hamiltonian systems is "SymplecticPartitionedRungeKutta". The option would look something like this:

Method -> {"SymplecticPartitionedRungeKutta", "DifferenceOrder" -> 10, "PositionVariables" -> {x[t]}} 

The one requirement is that the Hamiltonian be separable, that is, of the form $H(p,q,t)==T(p)+V(q,t)$. There is a bug, unfortunately still present in V10.4.1, for which a fix can be found in my answer here: SymplecticPartitionedRungeKutta shows strange error. Using the function cs[] found in it, we can invoke the "SymplecticPartitionedRungeKutta" method as follows. It turns out that if the default StartingStepSize is reduced to period/1000, we get excellent accuracy.

Block[{NDSolve`SPRKDump`CheckSeparability = cs}, Module[{m, x0, v0, u, f(* ... *)}, (* ...same as before... *) {ndsol} = NDSolve[sys, {x, p}, {t, 0, period}, Method -> {"SymplecticPartitionedRungeKutta", "DifferenceOrder" -> 10, "PositionVariables" -> {x[t]}}, InterpolationOrder -> All(*, StartingStepSize -> period/1000*)] ]] Plot[ham - ham0 /. t -> Mod[t, period] /. ndsol // RealExponent // Evaluate, {t, 0, period}, PlotPoints -> 201, PlotRange -> {-18, 0}, GridLines -> {None, {-7, -8}}]