Potential hill simulation

Question

I attempted to simulate the motion of a particle approaching to a potential hill. The particle should experience a force of -1/r^2, if the equation for the potential hill is 1/r. So the acceleration for the particle should be -c/r^2, where c is an arbitrary proportional constant. The particle is expected to bounce back when close enough to the hill, but it passed through the hill instead, despite of different values of c. Is there anything wrong with my code or the equations of motions. Secondly, when PlotPoints->10 changed to PlotPoints->30 in order to obtain an intact potential hill, the program becomes unbearable slow. How to get around this?

Clear[x, y, z, t, ux, uy, bx, by, cnt];
t = -2;
ux = 1; uy = 1;
bx = -1; by = -1;
obx = -1; oby = -1;
cnt = 0.2;
ax[x_, y_] := -cnt*x/( x^2 + y^2)^(3/2)  ;
ay[x_, y_] := -cnt*y/( x^2 + y^2)^(3/2)   ;

field = ContourPlot3D[
   1/(x^2 + y^2)^(1/2) - z == 0, {x, -2, 2}, {y, -2, 2}, {z, 0, 50}, 
   PlotPoints -> 10];

particle[t_] := (
   bx = ux t + 0.5 ax[obx, oby] t^2;
   by = uy t + 0.5 ay[obx, oby] t^2;
   bz = 1/(bx^2 + by^2)^(1/2) ;
   obx = bx; oby = by;
   Graphics3D[{
     Red, Sphere[{bx, by, bz}, 0.1]}]);

Dynamic[
 Refresh[
  Row[{Show[field, particle[t], ImageSize -> {300, 300}], 
    If[t < 2, t = t + 0.1, t = -2]}]
  , UpdateInterval -> 1
  ]
   ]

Answer 1

A couple things: You don't update velocity (obx and oby appear to be unused?). You do your timestep wrong (you want x+=velocity*timestep, not x=velocity*time!). You also have your direction of your force wrong. You should take the negative gradient of potential. The gradient of 1/r is -1/r^2. The negative of that is 1/r^2. So you get the common sense result - it pushes away from the center.

The usual way to do this, though, is to use NDSolve:

soln = First@NDSolve[{x''[t] == x[t]/(x[t]^2 + y[t]^2)^(3/2),
     y''[t] == y[t]/(x[t]^2 + y[t]^2)^(3/2),
     x[0] == -1,
     y[0] == -1,
     x'[0] == 1.8,
     y'[0] == 2}, {x, y}, {t, 0, 2}];

plt = Plot3D[1/Sqrt[x^2 + y^2], {x, -2, 2}, {y, -2, 2}, 
   PlotRange -> {-1, 5}];

Animate[Show[plt, 
   Graphics3D[{Red, 
     Sphere[{x[t], y[t], 1/Sqrt[x[t]^2 + y[t]^2]}, 0.1]}],BoxRatios->Automatic] /. 
  soln, {t, 0, 2}]