RST v1.1 – 1D Substrate–Resonance Simulator (Leapfrog + Damping)

January 18, 2026

RST v1.1 – 1D Substrate–Resonance Simulator (Leapfrog + Damping)

This Python script implements a 1D Reactive Substrate Theory (RST) simulator with a real substrate field S(x,t) and a complex resonance field Ψ(x,t). It uses a leapfrog-style time integrator, Gaussian initial seeds, a damping phase to let soliton-like structures relax, and a simple absorbing boundary layer to reduce reflections at the edges.


import numpy as np
import matplotlib.pyplot as plt

# -----------------------------
# PARAMETERS
# -----------------------------
# Field couplings
beta = 1.0       # S self-interaction
lambda_ = 1.0    # Psi self-interaction
mu = 1.0         # Psi mass term
kappa = 0.5      # S-Psi coupling

# Grid and time
dx = 0.05
dt = 0.01
xmax = 20.0
tmax = 10.0

# Damping for soliton relaxation
damping_steps = 300
damping_factor = 0.98

# Absorbing boundary width
absorb_width = int(0.1 * (xmax/dx))  # last 10% of grid

# -----------------------------
# GRID
# -----------------------------
x = np.arange(0, xmax, dx)
nx = len(x)

# Fields: S real, Psi complex
S = np.zeros(nx)
Psi = np.zeros(nx, dtype=complex)

# Velocities for leapfrog
S_dot = np.zeros(nx)
Psi_dot = np.zeros(nx, dtype=complex)

# -----------------------------
# INITIAL CONDITIONS (Gaussian Seeds)
# -----------------------------
x0 = xmax/2
S_amp, S_sigma = 1.0, 1.0
Psi_amp, Psi_sigma = 0.5, 1.0

S = S_amp * np.exp(-(x - x0)**2 / (2 * S_sigma**2))
Psi = Psi_amp * np.exp(-(x - x0)**2 / (2 * Psi_sigma**2))

# -----------------------------
# HELPER FUNCTIONS
# -----------------------------
def laplacian(f, dx):
    lap = np.zeros_like(f)
    lap[1:-1] = (f[2:] - 2*f[1:-1] + f[:-2]) / dx**2
    lap[0] = lap[1]
    lap[-1] = lap[-2]
    return lap

def apply_absorbing_boundary(vel, width):
    vel[:width] *= np.linspace(0,1,width)
    vel[-width:] *= np.linspace(1,0,width)
    return vel

def compute_energy(S, Psi, dx):
    kinetic_S = 0.5 * (np.gradient(S, dx))**2
    potential_S = 0.25 * beta * S**4
    kinetic_Psi = (np.abs(np.gradient(Psi, dx)))**2
    potential_Psi = 0.5 * mu * np.abs(Psi)**2 + 0.25 * lambda_ * np.abs(Psi)**4
    coupling = kappa * S * np.abs(Psi)**2
    energy_density = kinetic_S + potential_S + kinetic_Psi + potential_Psi + coupling
    return np.sum(energy_density) * dx

# -----------------------------
# TIME EVOLUTION (Leapfrog)
# -----------------------------
n_steps = int(tmax/dt)

for step in range(n_steps):
    # 1) Compute accelerations
    S_ddot = laplacian(S, dx) - beta*S**3 - kappa*np.abs(Psi)**2
    Psi_ddot = laplacian(Psi, dx) - mu*Psi - lambda_*np.abs(Psi)**2*Psi - kappa*S*Psi
    
    # 2) Half-step velocity update
    S_dot += 0.5 * dt * S_ddot
    Psi_dot += 0.5 * dt * Psi_ddot
    
    # 3) Full-step field update
    S += dt * S_dot
    Psi += dt * Psi_dot
    
    # 4) Recompute accelerations at new step
    S_ddot = laplacian(S, dx) - beta*S**3 - kappa*np.abs(Psi)**2
    Psi_ddot = laplacian(Psi, dx) - mu*Psi - lambda_*np.abs(Psi)**2*Psi - kappa*S*Psi
    
    # 5) Half-step velocity update
    S_dot += 0.5 * dt * S_ddot
    Psi_dot += 0.5 * dt * Psi_ddot
    
    # 6) Damping during soliton relaxation
    if step < damping_steps:
        S_dot *= damping_factor
        Psi_dot *= damping_factor
    
    # 7) Apply absorbing boundaries
    S_dot = apply_absorbing_boundary(S_dot, absorb_width)
    Psi_dot = apply_absorbing_boundary(Psi_dot, absorb_width)
    
    # 8) Diagnostics and plotting every 200 steps
    if step % 200 == 0:
        E_tot = compute_energy(S, Psi, dx)
        plt.figure(figsize=(8,3))
        plt.plot(x, S, label='S')
        plt.plot(x, np.abs(Psi), label='|Psi|')
        plt.title(f"Step {step}, Total Energy = {E_tot:.3f}")
        plt.legend()
        plt.show()

# -----------------------------
# FINAL ENERGY
# -----------------------------
E_tot = compute_energy(S, Psi, dx)
print("Final Total Energy:", E_tot)

You can paste this block directly into a Blogger post (in HTML mode). Readers familiar with Python and numerical PDEs can copy the code into a local file (e.g., rst_1d_sim.py), install numpy and matplotlib, and start experimenting with RST soliton-like dynamics.

Conspiranon