Modelado de un pozo de potencial

ACTUALIZADO: ver más abajo.

Sus entradas NDSolve parecen estar haciendo lo que esperaría para una masa alrededor de un centro gravitacional. Usando:

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[0] == 1, 
      y[0] == 0, x'[0] == 0, y'[0] == 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}}]

donde el centro gravitacional está en 0,0 y la masa está en (x,y)=(1,0) en t=0, obtengo una serie de órbitas limitadas y (probablemente) ilimitadas de varios tamaños en un rango de inicial velocidades:

Como dijo @alemi, no habrá una espiral a menos que incluya un término de amortiguación. Las órbitas para este tipo de potencial están libres (de entrada y luego de salida) o son bucles cerrados.

ACTUALIZACIÓN: la fuerza en las ecuaciones de movimiento tal como están escritas siempre es radial desde (x,y)=(0,0). La fuerza debe estar en la dirección (a,b)-(x,y). Las ecuaciones de movimiento correctas son:

{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[0] == -2, y[0] == -2, x'[0] == sp, y'[0] == sp}

Observe el cambio de 'a' y 'b' en los numeradores , no solo el cálculo de la distancia en el denominador. Las ecuaciones funcionaron en mi prueba porque coloqué el pozo de potencial en (0,0) y moví la masa fuera del eje.

Si tu potencial es$\propto 1/r$, estás simulando efectivamente la gravedad. Si te da una mejor intuición, imagínalo como la tierra alrededor del sol.

Siempre que su solucionador numérico esté haciendo un trabajo decente, no debe esperar que la bola entre en espiral, la solución correcta sería una sección cónica, es decir, orbitaría el origen en una trayectoria elíptica si las condiciones iniciales equivalen a un energía ligada, o hiperbólica si las condiciones iniciales equivalen a una energía no ligada. Los bocetos que publica sugieren que se encuentra en el caso no vinculado, por lo que, incluso para ver las órbitas, debe disminuir la magnitud de sus velocidades iniciales.

Por supuesto, esto solo es cierto siempre que esté usando un integrador simpléctico (uno que conserva energía, que imagino que Mathematica usa de manera predeterminada). Un integrador débil como el simple Euler tendría la bola en espiral, pero esto es puramente un efecto de la pobre integración.

Si realmente desea que entre en espiral, debe agregar algún tipo de término de amortiguación viscosa

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

o algo así como la resistencia del aire, una fuerza de la forma

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

Además, después de la actualización, está claro que todavía está disparando casi directamente al origen con una velocidad demasiado alta. Es como si estuvieras modelando dispersión gravitacional, no órbitas. Estás variando tu velocidad inicial , no tu velocidad inicial . Para ser completamente general, ¿por qué no establece

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

Y juegue con los controles deslizantes v0y thpor separado para tener una idea de lo que está sucediendo. Aquí la thvariable controlará la dirección inicial de su velocidad inicial y v0controlará la velocidad inicial.

En cuanto al valor que debe tener la velocidad inicial, puede considerar la energía inicial para tener una idea de la magnitud adecuada. En$t=0$, tienes

$$E = \frac 12 v^2 -\frac{1}{r}$$ Si esta energía es positiva, ignorando el amortiguamiento, la partícula no está unida y debería dispararse hasta el infinito. Intente usar la configuración inicial que Jason A usó en su respuesta. Ponga el centro del potencial en cero, haga que la partícula comience en $x(0) = 1, y(0) = 0$y establecer $\dot x(0) = 0, \dot y(0) = v$, dónde $v$está cerca de 1 más o menos.