# Solving equation of motion for position dependent force

Consider the functions

U[x_] = 1/x^2 - 2/x;
F[x_] = -1*U'[x];


But here comes the tricky part:

Write down the object's equation of motion, i.e., Newton's second law, in the form of a differential equation involving x[t] and its derivatives. To solve this equation numerically we use NDSolve. Complete the code for NDSolve below. Let the particle begin from rest at x[0] = x0 = 0.55. It has a mass of m = 0.01.

m = 0.01; x0 = 0.55; v0 = 0;
x[t_] := (F[x]/(2*m))*t^2 + x'[t]*t + x0
ClearAll[x];
solution = NDSolve[{x[t] == x'[t] + x0, x'[0] = 0 }, x, {t, 0, 18}];
x = x /. solution[[1]];

• So you define x[t_] := (F[x]/(2*m))*t^2 + x'[t]*t + x0 then you ClearAll[x]? This clears you definition, or am I lost? – mattiav27 Aug 7 '16 at 11:52
• Also if this is an homework, you should the appropriate tag... – mattiav27 Aug 7 '16 at 11:53
• Cool. Not homework. I have free time and am trying to learn Mathematica. Found this Physics problem but it seems a bit hard for me. Didn't realize ClearAll[x] does this. Just saw it being used... – Dylan Solms Aug 7 '16 at 12:03
• you obviously don't want to define x before you solve for it anyway. Your equation is second order F=m x''[t], hope that helps. (I assume F is force, though you don't actually say that ). I,ll add your attempt looks like you are trying to make some use of some constant acceleration expressions. Your acceleration is most definately not constant, you cant use those. – george2079 Aug 7 '16 at 12:56
• So, you are trying to solve a differential equation where a particle undergoes a position dependent force? – Feyre Aug 7 '16 at 12:59

Set the force:

u = x[t]^-2 - 2/x[t];
f = -D[u, x[t]];


Initial conditions:

{m,x0,v0} = {0.01,0.55,0};


Solve and plot:

s = NDSolve[{x''[t] == f/m, x[0] == x0, x'[0] == v0},
x[t], {t, 0, 3}];
Plot[x[t] /. s, {t, 0, 3}]


You can find the period with a FindMinimum[]

FindMinimum[x[t] /. s, {t, 2.3}]


{0.55, {t -> 2.33751}}

Note, that this is the same as the initial condition, the velocity is near to 0:

D[x[t] /. s, t] /. %[[2, 1]]


-9.42075*10^-6

As @michael-e2 suggested, the total energy is negative.

Plotting the total energy over three periods reveals a small numerical error, but the energy is clearly conserved:

Plot[(0.5 m D[x[t] /. s, t]^2 + u /. s) /. t -> ti, {ti, 0,
3 t /. %5[[2]]}, PlotRange -> {-0.33058, -0.330578},
Ticks -> {Automatic, {-0.330579, -0.330578}}]


And for the fun:

b = Table[
VectorPlot3D[{0, 0, 2/x^3 - 2/x^2}, {a, -1, 1}, {b, -1, 1}, {x, i,
6}, VectorPoints -> {3, 3, 9}], {i, -0.5, 1, 1.5}];
Manipulate[
Show[{Graphics3D[{Sphere[{0, 0, Evaluate[x[t] /. s][[1]]}, 0.4]},
PlotRange -> {{-1, 1}, {-1, 1}, {-0, 6}}] /. t -> tl}[[1]],
b], {tl, 0, 2.3, 0.05}]


• What is the meaning of arrows on the last plot? Is it possible to make them 3d to be consistent with the 3d ball? – yarchik Aug 7 '16 at 16:01
• @yarchik I added a VectorPlot3D of the force, those are the default arrows. I could, but I threw away the code and it's dinner time. Feel free to change the gif if you feel like it. – Feyre Aug 7 '16 at 17:23
• @Feyre I think you should include the code for the fun part. Pretty please. – dearN Aug 7 '16 at 21:15
• @drN Because you asked so nicely. – Feyre Aug 8 '16 at 8:39

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).

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[{NDSolveSPRKDumpCheckSeparability = 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}}]


• Nice summary of what works and doesn't work in the Hamiltonian formulation - probably overkill for this question, though... I felt tempted to close it, but now I can't bring myself to vote for that. – Jens Aug 8 '16 at 3:23
• I learned a lot. +1 – LLlAMnYP Aug 8 '16 at 8:53
• "x[t] and its derivatives" I take that to mean standard classical mechanics. Still, such a long post makes me feel rather outdone. – Feyre Aug 8 '16 at 11:51
• @Feyre, Jens, et al. -- Feyre's original bouncing ball is what inspired me. And I've been wanting to work up a Hamiltonian mechanics example for class, and this seemed a like a good choice. It's like a Lennard-Jones potential, which I might use instead. Why not share my explorations, even if they are overkill? (Sometimes, I don't have the time.) Also, the OP's "I have free time and am trying to learn Mathematica" was an inspiration, too, more so than "x[t]` and its derivatives", obviously. :) – Michael E2 Aug 8 '16 at 16:13
• @MichaelE2 I agree, we have to be able to enjoy coding here. – Feyre Aug 8 '16 at 17:27