Central Force Problems

da PowerShell -- jupyter nbconvert --to html notebook.ipynb

Moto radiale

L’esercizio simula numericamente il moto radiale di una particella soggetta a una forza centrale gravitazionale newtoniana.

Si considera un potenziale gravitazionale $$\large V(r) = -\frac{1}{r}, $$ e si costruisce il potenziale efficace $$\large V_{\text{eff}}(r) = V(r) + \frac{L^2}{2mr^2}, $$ che tiene conto del momento angolare conservato $L=m\, r^2\, \dot{\theta}$.

L’equazione del moto radiale è $$\large m \ddot r = -\frac{dV_{\text{eff}}}{dr}, $$ dove la derivata del potenziale efficace viene calcolata numericamente mediante differenze finite centrate.

Il sistema di equazioni del primo ordine $$\large \begin{cases} \dot r = v \\ \dot v = -\frac{1}{m}\frac{dV_{\text{eff}}}{dr} \end{cases} $$ viene integrato nel tempo tramite il metodo di Runge--Kutta del quarto ordine (RK4), a partire da condizioni iniziali assegnate per $r(0)$ e $\dot r(0)$.

Infine, viene rappresentata graficamente l’evoluzione temporale della coordinata radiale $r(t)$, evidenziando il comportamento orbitale determinato dal potenziale efficace.

import numpy as np
import matplotlib.pyplot as plt 

def gravitational_potential(r):
    """
    Gravitational potential energy function V(r) = -G * ml * m2 / r.
        :param r: Radial distance.
    :return: 
        Gravitational potential energy.
    """
    return -1.0 / r

def effective_potential(r, L, m, V_func):
    """
    Calculate the effective potential.
        :param r: Radial distance.
        :param L: Angular momentum.
        :param m: Mass of the particle.
        :param V_func: Potential energy function V(r).
    :retum: 
        Effective potential value.
    """
    return V_func(r) + L**2 / (2 * m * r**2)

def dV_eff_dr(r, L, m, V_func):
    """
    Calculate the derivative of effective potential with respect to r.
        :param r: Radial distance.
        :param L: Angular momentum.
        :param m: Mass of the particle.
        :param V_func: Potential energy function V(r).
    :retum: 
        Derivative of the effective potential.
    """
    eps = 1e-6
    return (
        effective_potential(r + eps, L, m, V_func)
        - effective_potential(r - eps, L, m, V_func)
    ) / (2 * eps)

def radial_acceleration(r, L, m, V_func):
    return - dV_eff_dr(r, L, m, V_func) / m

def central_force_solver_RK4(r0, dr0, m, L, V_func, t_max, dt):
    """
    Numerical integration of radial motion using Euler's method.
        :param r0: Initial radial position.
        :param dr0: Initial radial velocity.
        :param m: Mass of the particle.
        :param L: Angular momentum.
        :param V_func: Potential energy function V(r).
        :param t_max: Total time for simulation.
        :param dt: Time step.
    :return: 
        Arrays of time, radial positions, and velocities.
    """
    def acc(r):
        return radial_acceleration(r, L, m, V_func)

    t = 0.0
    r = r0
    v = dr0

    time = []
    positions = []
    velocities = []

    while t < t_max:

        # RK4
        k1_r = v
        k1_v = acc(r)

        k2_r = v + 0.5 * dt * k1_v
        k2_v = acc(r + 0.5 * dt * k1_r)

        k3_r = v + 0.5 * dt * k2_v
        k3_v = acc(r + 0.5 * dt * k2_r)

        k4_r = v + dt * k3_v
        k4_v = acc(r + dt * k3_r)

        r += dt * (k1_r + 2*k2_r + 2*k3_r + k4_r) / 6
        v += dt * (k1_v + 2*k2_v + 2*k3_v + k4_v) / 6

        time.append(t)
        positions.append(r)
        velocities.append(v)

        t += dt

    return np.array(time), np.array(positions), np.array(velocities)

r0    = 2.5
dr0   = -0.3
m     = 1.0
L     = 0.7 # momento angolare = m*r^2*rheta_punto
t_max = 20.0
dt = 0.0002


time, r, v = central_force_solver_RK4(
    r0, dr0, m, L, gravitational_potential, t_max, dt
)

plt.figure(figsize=(12, 3))
plt.plot(time, r)
plt.xlabel("t")
plt.ylabel("r(t)")
plt.grid()
plt.show()

E = 0.5 * m * v**2 + effective_potential(r, L, m, gravitational_potential)

plt.figure(figsize=(12, 3))
plt.ylim([-1, 1])
plt.plot(time, E-E.max())
plt.xlabel("t")
plt.ylabel("Variazione energia totale")
plt.grid()
plt.show()