Modellierung eines Potentialbrunnens

AKTUALISIERT: Siehe unten.

Ihre NDSolve-Eingaben scheinen das zu tun, was ich für eine Masse um ein Gravitationszentrum erwarten würde. Verwendung:

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}}]

wo das Gravitationszentrum bei 0,0 und die Masse bei (x,y)=(1,0) bei t=0 ist, erhalte ich eine Reihe von begrenzten und (wahrscheinlich) unbegrenzten Umlaufbahnen unterschiedlicher Größe über einen Anfangsbereich Geschwindigkeiten:

Wie @alemi feststellte, wird es keine Spirale geben, es sei denn, Sie fügen einen Dämpfungsbegriff ein. Die Umlaufbahnen für diese Art von Potenzial sind entweder ungebunden (Einflug, dann Ausflug) oder geschlossene Schleifen.

UPDATE: Die Kraft in den geschriebenen Bewegungsgleichungen ist immer radial von (x,y)=(0,0). Die Kraft sollte in Richtung (a,b)-(x,y) wirken. Die richtigen Bewegungsgleichungen lauten:

{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}

Beachten Sie die Verschiebung um 'a' und 'b' in den Zählern , nicht nur die Entfernungsberechnung im Nenner. Die Gleichungen funktionierten in meinem Test, weil ich den Potentialtopf bei (0,0) platzierte und die Masse von der Achse wegbewegte.

Wenn Ihr Potenzial ist$\propto 1/r$, simulieren Sie effektiv die Schwerkraft. Wenn es Ihnen eine bessere Intuition gibt, stellen Sie es sich als die Erde um die Sonne vor.

Solange Ihr numerischer Löser gute Arbeit leistet, sollten Sie nicht erwarten, dass sich die Kugel spiralförmig hineindreht, die richtige Lösung wäre ein Kegelschnitt, dh sie würde den Ursprung auf einer Ellipsenbahn umkreisen, wenn die Anfangsbedingungen gleich a sind gebundene Energie oder hyperbolisch, wenn die Anfangsbedingungen einer ungebundenen Energie entsprechen. Die von Ihnen geposteten Skizzen deuten darauf hin, dass Sie sich im ungebundenen Fall befinden. Um also Umlaufbahnen zu sehen, müssen Sie die Größe Ihrer Anfangsgeschwindigkeiten verringern.

Dies gilt natürlich nur, solange Sie einen symplektischen Integrator verwenden (einen, der Energie spart, den Mathematica meiner Meinung nach standardmäßig verwendet). Ein schwacher Integrator wie ein einfacher Euler hätte die Kugelspirale drin, aber das ist nur ein Effekt der schlechten Integration.

Wenn Sie wirklich wollen, dass es spiralförmig eindringt, müssen Sie eine Art viskosen Dämpfungsterm hinzufügen

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

oder so etwas wie Luftwiderstand, eine Kraft der Form

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

Außerdem ist nach dem Update klar, dass Sie immer noch mit zu hoher Geschwindigkeit fast direkt auf den Ursprung schießen. Es ist, als würden Sie Gravitationsstreuung modellieren, nicht Umlaufbahnen. Sie variieren Ihre Anfangsgeschwindigkeit , nicht Ihre Anfangsgeschwindigkeit . Um ganz allgemein zu sein, warum setzen Sie nicht

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

Spielen Sie separat mit den Schiebereglern v0und herum th, um ein Gefühl dafür zu bekommen, was vor sich geht. Hier thsteuert die Variable die Anfangsrichtung Ihrer Anfangsgeschwindigkeit und v0die Anfangsgeschwindigkeit.

Was den Wert der Anfangsgeschwindigkeit betrifft, können Sie die Anfangsenergie berücksichtigen, um ein Gefühl für die entsprechende Größe zu bekommen. Bei$t=0$, du hast

$$E = \frac 12 v^2 -\frac{1}{r}$$ Wenn diese Energie positiv ist und die Dämpfung außer Acht gelassen wird, ist das Teilchen ungebunden und sollte einfach ins Unendliche schießen. Versuchen Sie, die Anfangseinstellungen zu verwenden, die Jason A in seiner Antwort verwendet hat. Setzen Sie den Mittelpunkt des Potentials auf Null, lassen Sie das Teilchen bei beginnen $x(0) = 1, y(0) = 0$und einstellen $\dot x(0) = 0, \dot y(0) = v$, Wo $v$ist nahe 1 oder so.