Αλγόριθμος Velocity Verlet

import numpy as np
import matplotlib.pyplot as plt

# ==============================
# ΠΑΡΑΜΕΤΡΟΙ ΠΡΟΣΟΜΟΙΩΣΗΣ
# ==============================

dt = 0.01      # χρονικό βήμα Δt 
steps = 200    # αριθμός βημάτων, μπορεί να αλλάξει για μικρότερη ή μεγαλύτερη εξομοίωση
m = 1.0        # μάζα σωματιδίων

# Παράμετροι δυναμικού Lennard–Jones
epsilon = 1.0  # βάθος δυναμικού 
sigma = 1.0    # χαρακτηριστική απόσταση (όπου V=0)

# ==============================
# ΑΡΧΙΚΕΣ ΣΥΝΘΗΚΕΣ
# ==============================

# αρχικές Θέσεις δύο σωματιδίων (σε 2D)
r = np.array([
    [-1.5, 0.0],
    [ 1.5, 0.0]
])

# αρχικές ταχύτητες δύο σωματιδίων
v = np.array([
    [0.5, 0.2],
    [-0.5, -0.2]
])

# ==============================
# ΔΥΝΑΜΗ LENNARD-JONES
# ==============================

def lennard_jones_forces(r):
    """
    Υπολογίζει τις δυνάμεις μεταξύ όλων των ζευγών σωματιδίων.
    
    Επιστρέφει πίνακα forces ίδιου σχήματος με το r.
    """
    N = len(r)                    # αριθμός σωματιδίων
    forces = np.zeros_like(r)     # αρχικοποίηση δυνάμεων

    # Υπολογίζουμε μόνο για i < j (κάθε ζεύγος μία φορά)
    for i in range(N):
        for j in range(i+1, N):

            # διάνυσμα απόστασης r_i - r_j
            rij = r[i] - r[j]

            # μέτρο απόστασης |r| σωματιδίων
            dist = np.linalg.norm(rij)

            # αποφυγή διαίρεσης με το μηδέν
            if dist == 0:
                continue

            # υπολογισμός όρων (σ/r)^6 και (σ/r)^12
            sr6 = (sigma / dist)**6
            sr12 = sr6**2

            # μέτρο δύναμης (scalar)
            force_mag = 24 * epsilon * (2*sr12 - sr6) / dist**2

            # διάνυσμα δύναμης
            fij = force_mag * rij

            # δράση-αντίδραση 
            forces[i] += fij
            forces[j] -= fij

    return forces

# ==============================
# ΑΡΧΙΚΗ ΕΠΙΤΑΧΥΝΣΗ
# ==============================

# a = F/m
a = lennard_jones_forces(r) / m

# ==============================
# ΑΠΟΘΗΚΕΥΣΗ ΤΡΟΧΙΑΣ
# ==============================

# λίστα που κρατά όλες τις προηγούμενες θέσεις (για τα plots)
traj = [[] for _ in range(len(r))]

# ==============================
# ΕΝΕΡΓΟΠΟΙΗΣΗ INTERACTIVE PLOT
# ==============================

plt.ion()  # επιτρέπει live animation

# ==============================
# ΚΥΡΙΟ LOOP ΠΡΟΣΟΜΟΙΩΣΗΣ
# ==============================

for step in range(steps):

    # ---------------------------------
    # ΒΗΜΑ 1: μισό βήμα στην ταχύτητα
    # v(t + dt/2) = v(t) + 1/2 a(t) dt
    # ---------------------------------
    v_half = v + 0.5 * a * dt

    # ---------------------------------
    # ΒΗΜΑ 2: ενημέρωση θέσης
    # r(t + dt) = r(t) + v_half * dt
    # ---------------------------------
    r = r + v_half * dt

    # ---------------------------------
    # ΒΗΜΑ 3: νέες δυνάμεις -> νέα επιτάχυνση
    # ---------------------------------
    a_new = lennard_jones_forces(r) / m

    # ---------------------------------
    # ΒΗΜΑ 4: ολοκλήρωση ταχύτητας
    # v(t + dt) = v_half + 1/2 a(t+dt) dt
    # ---------------------------------
    v = v_half + 0.5 * a_new * dt

    # ενημέρωση επιτάχυνσης για το επόμενο βήμα
    a = a_new

    # ---------------------------------
    # ΑΠΟΘΗΚΕΥΣΗ ΤΡΟΧΙΑΣ
    # ---------------------------------
    for i in range(len(r)):
        traj[i].append(r[i].copy())

    # ==============================
    # PLOT
    # ==============================

    plt.clf()  # καθαρίζει το προηγούμενο frame

    # σχεδίαση τροχιών 
    for i in range(len(r)):
        t = np.array(traj[i])
        plt.plot(t[:,0], t[:,1], label=f"particle {i}")

    # σχεδίαση σωματιδίων 
    for i in range(len(r)):
        circle = plt.Circle((r[i,0], r[i,1]), 0.2)
        plt.gca().add_patch(circle)

    # όρια άξονων
    plt.xlim(-5, 5)
    plt.ylim(-5, 5)

    # σωστή αναλογία αξόνων (ίδια κλίμακα x και y)
    plt.gca().set_aspect('equal')

    plt.title(f"Step {step}")
    plt.grid()

    # μικρή παύση για το animation
    plt.pause(0.01)


plt.ioff()
plt.show()