Implementierung eines Kugelpendels

Das Problem scheint zu sein, dass Sie die Schrittgröße beim Integrieren nicht rechtzeitig berücksichtigen. Dies sollte offensichtlich sein, wenn Sie Dinge wie tun

phi_v += phi_f

im Code. ¹ Die Beschleunigung kann nicht einfach zur Geschwindigkeit addiert werden (die Einheiten stimmen nicht überein!). Die erwartete Beziehung ist,

$$a_\phi =\frac{\mathrm d v_\phi}{\mathrm dt}\approx\frac{\Delta v_\phi}{\Delta t}\implies v_\phi\approx v_\phi'+ a_\phi \Delta t$$ Wo$v_\phi'$ist der vorherige Wert ist die Geschwindigkeit.

Was Sie tun sollten, ist die Verwendung der 4 Gleichungen (2 Positionen, 2 Geschwindigkeiten),

$$\mathbf{v}=\frac{\mathrm d\mathbf x}{\mathrm dt} \text{ and } \mathbf{a}=\frac{1}{m}\frac{\mathrm d\mathbf f}{\mathrm dt} ,$$ Wo$\mathbf v=[\dot{\phi},\,\dot{\theta}]^T$und ähnlich für$\mathbf x$, und lösen Sie diese über einen Leap-Frog-Integrator wie Velocity Verlet . Da Sie jedoch eine geschwindigkeitsabhängige Beschleunigung haben (dh$a= f(\mathbf x,\,\mathbf v)$), müssen Sie eine modifizierte Version des Integrators verwenden (worüber ich in dieser anderen Antwort von mir spreche ),

$$\begin{align} a_1(t)&=\text{compute from $x(t)$ and $v(t)$}\tag{1}\\ x(t+\Delta t)&=x(t)+v(t)\Delta t+\frac{1}{2}a_1(t)\Delta t^2\tag{2}\\ v(t+\Delta t)&=v(t)+a_1(t)\Delta t\tag{3}\\ a_2(t+\Delta t)&=\text{compute from $x(t+\Delta t)$, $v(t+\Delta t)$}\tag{4}\\ v(t+\Delta t)&=v(t+\Delta t)+\frac{1}{2}\left(a_2(t+\Delta t)-a_1(t)\right)\Delta t\tag{5} \end{align}$$ die Sie dann einfach durchlaufen, bis Sie zufrieden sind (z. B.$t\geq t_\text{end}$). Bereitgestellt$\Delta t$klein genug ist, sollte dies einen stabileren Löser liefern als die Geschwindigkeits-Verlet-Integration allein.

^{1. OP hat den timestepsTerm zu den Gleichungen in Version 6 hinzugefügt ; Diese Antwort wurde vor dieser Bearbeitung gepostet.}

Es sieht für mich so aus, als ob Sie versuchen, ein System gewöhnlicher Differentialgleichungen zu lösen, indem Sie die Lösung für kleine Schritte schätzen. Hab ich recht? Wenn ja, dann gibt es einiges zu beachten.

Erstens benötigen Sie eine gute Schätzung für den Ableitungsoperator, um den zukünftigen Zustand unter den Anfangsbedingungen zu approximieren. Diese Art von Problem wird als Anfangswertproblem (ivp) bezeichnet.

Da die Gleichung sowohl Zeitableitungen erster als auch zweiter Ordnung beinhaltet, besteht ein gemeinsamer Schritt darin, Impulsvariablen zu definieren und ein größeres System zu propagieren. Dies hilft, Fehler zu vermeiden, wenn die Gesamtzeit zunimmt. Ich habe dieses Problem bei Raytracing-Codes gesehen, und was ich vorschlage, ist so häufig, dass die meisten Leute niemals eine Gleichung zweiter Ordnung propagieren. Die Idee ist:

p_theta = d(theta)/dt

p_phi = d(phi)/dt

dann ist deine gleichung

d(p_theta)/dt = (p_phi)^2*cos(theta)*sin(theta) - g/l*sin(theta)

d(p_phi)/dt = -2*p_theta*p_phi*cot(theta)

plus die beiden Gleichungen, die die p definieren.

Ein einfacher Euler-Schritt würde implementiert als

p_theta(t0+dt) = p_theta(t0) + (d(p_theta)/dt)(t0)*dt
p_phi(t0+dt) = p_phi(t0) + (d(p_phi)/dt)(t0)*dt

plus die Gleichungen für Theta und Phi in den Definitionen von p.

Dies ist die typische Konfiguration. Nun, der Euler-Schritt ist sehr schlecht und wird niemals empfohlen. Sie sollten besser eine Methode höherer Ordnung wie RK4 oder RK5 (4) usw. mit Schrittgrößensteuerung implementieren.

Abgesehen davon ist die Verwendung von Winkeln manchmal ein Problem, da dies dazu führen kann, dass Sie bei einem Schritt den nächsten Wert nicht eindeutig bestimmen können, da das System singulär ist. In Luft- und Raumfahrtsimulationen verwenden sie Viertelionen, um dies zu beheben. Einzelheiten zur Mathematik oder Zipfel-Modellierung und -Simulation der Dynamik von Luft- und Raumfahrtfahrzeugen finden Sie in Goldstein Classical Mechanics. Ich denke, Sie brauchen diese Maschinerie jetzt nicht. Sie müssen die Gleichungen richtig schreiben und einen einfachen Schrittalgorithmus ausprobieren, bevor Sie zu anspruchsvoll werden. Ich würde denken, dass Python ein ODE-Solver-Paket hat, sodass Sie kein eigenes schreiben müssen, sondern nur richtig einrichten und aufrufen. Ich hoffe das hilft.