# Modeling a potential well

I attempted to simulate the interaction of a moving particle and a potential well in Mathematica. The particle should experience a force of -$1/r^2$, if the equation for the potential well is -$1/r$. The inputs for the NDSolve function are:

  x''[t] == -x[t]/((x[t] - a)^2 + (y[t] - b)^2)^(3/2),
y''[t] == -y[t]/((x[t] - a)^2 + (y[t] - b)^2)^(3/2),
x == -2,
y == -2,
x' == 2 ,
y' == 2,


A kind of spiraling motion or trapping should happen when the particle approaching to the well in a certain angle. But it does not happen in this simple simulation. The particle must experiences some centrifugal force when it is close to the well so that its direction is changed. What other terms should be added to represent this force in the NDSolve function? Any idea or any help would be much appreciated.

What I got is: What expected to see is: By adopting the suggestions given below, the terms for the NDSolve were now modified as follows:

soln[a_, b_, alp_, sp_] := Return[First@NDSolve[{
x''[t] == x[t]/((x[t] - a)^2 + (y[t] - b)^2)^(3/2) - alp*x'[t],
y''[t] == y[t]/((x[t] - a)^2 + (y[t] - b)^2)^(3/2) - alp*y'[t],
x == -2,
y == -2,
x' == sp,
y' == sp},
{x, y}, {t, 0, 3}]];


and sliders were added for adjusting. Still, no conic trajectory was observed. The particle moves as if in a straight line from the top view. • With those inputs you launch the particle directly into the divergence. I would expect then, due to numerical errors, it to overshoot (0, 0) slightly and come out the other way in a straight line. This is to say that the trajectory looks plausible given the inputs. If I'm right, were you to plot speed over time, you should get a non-constant value, and if you started the particle out at a slight angle, y' = 1, say, you should get the behavior you want (note that the starting speed might be too high for the particle to "notice" gravity). – alarge Aug 8 '14 at 2:32
• Actually, the position of the potential well is adjustable in the simulation, different values of initial speed and positions of the potential well also had been tried. No trace of any spiral had been observed. – user16069 Aug 8 '14 at 2:55
• In that case, this does not sound like a physics problem, but more like an issue with the particular software. There is a Mathematica.SE, and the folks there would probably be able to help better. – alarge Aug 8 '14 at 3:12
• @alarge His particular problem in this case is a physics one, he is simulating newtonian gravity and expecting it to spiral in. – alemi Aug 8 '14 at 3:15
• @alemi Agreed, there is a slight misunderstanding as to what the orbit should look like in the question as posed, but the main problem here seems to be that the software does not solve the equations: The trajectory is not deflected (expecting a hyperbola, a conic section, is reasonable but apparently does not happen despite varying the initial conditions). – alarge Aug 8 '14 at 3:18

UPDATED: See below.

Your NDSolve inputs seem to be doing what I would expect for a mass around a gravitational center. Using:

a = 0;
b = 0;
traj = Table[
s = NDSolve[{x''[t] == -x[t]/((x[t] - a)^2 + (y[t] - b)^2)^(3/2),
y''[t] == -y[t]/((x[t] - a)^2 + (y[t] - b)^2)^(3/2), x == 1,
y == 0, x' == 0, y' == v}, {x, y}, {t, -20, 20}];
{x[t], y[t]} /. s,
{v, 0.1, 2.1, 0.2}];
ParametricPlot[traj, {t, -20, 20}, PlotRange -> {{-3, 2}, {-4, 4}}]


where the gravitational center is at 0,0 and the mass is at (x,y)=(1,0) at t=0, I get a series of bounded and (probably) unbounded orbits of various sizes over a range of initial velocities: As @alemi stated, there will be no spiraling in unless you include a damping term. The orbits for this type of potential are either unbound (fly-in then fly-out) or are closed loops.

UPDATE: the force in the equations of motion as written is always radial from (x,y)=(0,0). The force should be in the direction (a,b)-(x,y). The correct equations of motion are:

{x''[t] == (a - x[t])/((x[t] - a)^2 + (b - y[t])^2)^(3/2) - alp*x'[t],
y''[t] == (b - y[t])/((x[t] - a)^2 + (y[t] - b)^2)^(3/2) -
alp*y'[t], x == -2, y == -2, x' == sp, y' == sp}


Notice the shift by 'a' and 'b' in the numerators, not just the distance calculation in the denominator. The equations worked in my test because I placed the potential well at (0,0) and moved the mass off axis.

• you can remove the (probably). $E = -1/r + \frac 12 v^2$, so in your case $v=\sqrt 2$ is the turning point. – alemi Aug 8 '14 at 3:50

If your potential is $\propto 1/r$, you're effectively simulating gravity. If it gives you better intuition, imagine it as the earth around the sun.

As long as your numerical solver is doing a decent job, you shouldn't expect the ball to spiral in, the correct solution would be a conic section, that is it would orbit the origin in an elliptical path if the initial conditions equate to a bound energy, or a hyperbolic one if the initial conditions equate to an unbound energy. The sketches you post suggest you are in the unbound case, so, even to see orbits you need to decrease the magnitude of your initial velocities.

This is of course only true long as you were using a symplectic integrator (one that conserves energy, which I imagine Mathematica uses by default). A weak integrator like simple Euler would have the ball spiral in, but this is a purely an effect of the poor integration.

If you truly want it to spiral in, you need to add some kind of viscous damping term

$$\vec F = -\alpha \vec v$$

or something like air resistance, a force of the form

$$\vec F = -\alpha v^2 \hat v= -\alpha |v| \vec v$$

Additionally, after you're update, it is clear that you are still shooting nearly straight at the origin with too high a speed. It is as if you are modelling gravitational scattering, not orbits. You are varying your initial speed, not your initial velocity. To be completely general, why don't you set

x' == v0 Cos[th]
y' == v0 Sin[th]


And play around with the v0 and th sliders separately to get a feel for what's going on. Here the th variable will control the initial direction of your initial velocity and v0 will control the initial speed.

As for what value the initial speed should be, you can consider the initial energy to get a sense of the appropriate magnitude. At $t=0$, you have

$$E = \frac 12 v^2 -\frac{1}{r}$$ If this energy is positive, ignoring the damping, the particle is unbound and should just shoot off to infinity. Try using the initial settings Jason A used in his answer. Put the center of the potential at zero, have the particle start at $x(0) = 1, y(0) = 0$ and set $\dot x(0) = 0, \dot y(0) = v$, where $v$ is near 1 or so.