Are the courses on CoursesHub.World free?

Yes, all courses on CoursesHub.World are completely free and open access. We believe in democratizing education and making university-level science courses available to everyone worldwide.

What subjects are covered on CoursesHub.World?

CoursesHub.World offers courses in physics (Quantum Mechanics, General Relativity, QFT, Plasma Physics, Cosmology), biology (Molecular Biology, Cell Physiology, Pharmacology), and earth sciences (Oceanography, Atmospheric Science, Climatology).

What level are the courses?

Our courses are designed at the graduate and advanced undergraduate level. They include rigorous mathematical derivations and are suitable for physics students, researchers, and serious self-learners.

Where do the video lectures come from?

Our 400+ video lectures come from world-renowned sources including MIT OpenCourseWare, Stanford, lectures by Nobel laureates, and other leading universities and educators.

← Back to Part V: Advanced Topics

Chapter 27: Numerical Relativity

Numerical relativity solves Einstein's equations computationally for situations too complex for analytical solutions. Essential for LIGO waveform templates and understanding binary black hole mergers.

BSSN Formulation

The Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation is the standard for stable evolution:

Key Variables

• Conformal factor: φ = ln(h)/12
• Conformal metric: γ̃_ij = e^-4φh_ij
• Traceless extrinsic curvature: Ã_ij
• Conformal connection functions: Γ̃ⁱ

Moving Punctures

The 2005 breakthrough: allow black hole "punctures" to move through the grid using "1+log" slicing and "gamma-driver" shift conditions.

1+log Slicing

∂_tα = -2αK (avoids singularity)

Gamma-Driver

∂_tβⁱ = 3/4 Bⁱ (coordinate motion)

Python: Wave Equation Toy Model

This simulation demonstrates finite difference methods used in numerical relativity by solving the 1D wave equation. The wave equation shares the same hyperbolic character as Einstein's equations.

Wave Equation Evolution Simulation

Python

Evolves a Gaussian pulse using finite differences with CFL stability

script.py172 lines

#!/usr/bin/env python3
"""
Numerical Relativity - 1D Wave Equation Evolution
Demonstrates finite difference methods used in numerical relativity
"""
import numpy as np
import matplotlib.pyplot as plt
import base64
from io import BytesIO

print("╔══════════════════════════════════════════════════════════════════╗")
print("║         NUMERICAL RELATIVITY: WAVE EQUATION TOY MODEL           ║")
print("╠══════════════════════════════════════════════════════════════════╣")
print("║  The 1D wave equation serves as a toy model for understanding   ║")
print("║  the hyperbolic nature of Einstein's equations                  ║")
print("╚══════════════════════════════════════════════════════════════════╝")

print("\n" + "="*70)
print("WAVE EQUATION AND DISCRETIZATION")
print("="*70)
print("\nWave equation: $ \\frac{\\partial^2 u}{\\partial t^2} = c^2 \\frac{\\partial^2 u}{\\partial x^2} $")
print("\nFinite difference approximation:")
print("$ \\frac{u_i^{n+1} - 2u_i^n + u_i^{n-1}}{\\Delta t^2} = c^2 \\frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\\Delta x^2} $")
print("\nRearranging for the update formula:")
print("$ u_i^{n+1} = 2u_i^n - u_i^{n-1} + \\left(\\frac{c\\Delta t}{\\Delta x}\\right)^2 (u_{i+1}^n - 2u_i^n + u_{i-1}^n) $")

# Parameters
nx = 200          # Number of spatial points
nt = 500          # Number of time steps
c = 1.0           # Wave speed
L = 1.0           # Domain length
dx = L / nx       # Spatial step
dt = 0.5 * dx / c # Time step (CFL condition: dt <= dx/c)

print("\n" + "="*70)
print("SIMULATION PARAMETERS")
print("="*70)
print(f"\n• Grid points (nx): {nx}")
print(f"• Time steps (nt): {nt}")
print(f"• Wave speed (c): {c}")
print(f"• Spatial step (Δx): {dx:.6f}")
print(f"• Time step (Δt): {dt:.6f}")
print(f"\nCFL number: $ \\sigma = \\frac{c \\Delta t}{\\Delta x} = {c*dt/dx:.2f} $")
print("\nCFL stability condition: $ \\sigma \\leq 1 $ ✓")

# Initial conditions: Gaussian pulse
x = np.linspace(0, L, nx)
sigma_gauss = 0.05
u = np.exp(-((x - 0.5)**2) / (2 * sigma_gauss**2))
u_prev = u.copy()  # Zero initial velocity

print("\n" + "="*70)
print("INITIAL CONDITIONS")
print("="*70)
print("\nGaussian initial data:")
print("$ u(x,0) = \\exp\\left(-\\frac{(x-0.5)^2}{2\\sigma^2}\\right) $")
print(f"\nGaussian width σ = {sigma_gauss}")
print("Initial velocity: $ \\frac{\\partial u}{\\partial t}\\bigg|_{t=0} = 0 $")

# Storage for snapshots
snapshots = []
snapshot_times = [0, 100, 200, 300, 400]

# Evolution
print("\n" + "="*70)
print("TIME EVOLUTION")
print("="*70)
print("\nEvolving wave equation...")

for n in range(nt + 1):
    if n in snapshot_times:
        snapshots.append((n * dt, u.copy()))
        print(f"  Snapshot at t = {n*dt:.3f} (step {n})")

if n < nt:
        # Compute new timestep
        u_new = np.zeros(nx)
        # Interior points
        courant = (c * dt / dx)**2
        for i in range(1, nx - 1):
            u_new[i] = 2*u[i] - u_prev[i] + courant * (u[i+1] - 2*u[i] + u[i-1])
        # Boundary conditions (fixed ends)
        u_new[0] = 0
        u_new[-1] = 0
        # Update
        u_prev = u.copy()
        u = u_new.copy()

print("\nEvolution complete!")

# Create visualization
fig, axes = plt.subplots(2, 3, figsize=(14, 8))
axes = axes.flatten()

colors = ['#4ade80', '#22d3ee', '#f472b6', '#facc15', '#a78bfa']

for idx, (t, u_snap) in enumerate(snapshots):
    ax = axes[idx]
    ax.plot(x, u_snap, color=colors[idx], linewidth=2)
    ax.fill_between(x, u_snap, alpha=0.3, color=colors[idx])
    ax.set_xlabel('x', fontsize=11)
    ax.set_ylabel('u(x,t)', fontsize=11)
    ax.set_title(f't = {t:.3f}', fontsize=12, fontweight='bold', color=colors[idx])
    ax.set_xlim(0, 1)
    ax.set_ylim(-1.2, 1.2)
    ax.grid(True, alpha=0.3)
    ax.axhline(y=0, color='gray', linestyle='-', alpha=0.5)

# Energy conservation check
ax = axes[5]
energies = []
for t, u_snap in snapshots:
    # Approximate energy: E ~ ∫(u_t² + c²u_x²)dx
    # For simplicity, just compute ∫u²dx
    E = np.trapz(u_snap**2, x)
    energies.append(E)

times_snap = [t for t, _ in snapshots]
ax.bar(range(len(energies)), energies, color=colors, alpha=0.7)
ax.set_xlabel('Snapshot', fontsize=11)
ax.set_ylabel('∫u²dx', fontsize=11)
ax.set_title('Energy Conservation Check', fontsize=12, fontweight='bold')
ax.set_xticks(range(len(energies)))
ax.set_xticklabels([f't={t:.2f}' for t in times_snap], fontsize=9)
ax.grid(True, alpha=0.3, axis='y')

plt.suptitle('Wave Equation Evolution (Toy Model for Numerical Relativity)',
             fontsize=14, fontweight='bold', color='white')
fig.patch.set_facecolor('#1a1a2e')
for ax in axes:
    ax.set_facecolor('#16213e')
    ax.tick_params(colors='white')
    ax.xaxis.label.set_color('white')
    ax.yaxis.label.set_color('white')
    ax.title.set_color('white')
    for spine in ax.spines.values():
        spine.set_color('gray')

plt.tight_layout()

# Save to base64
buf = BytesIO()
plt.savefig(buf, format='png', dpi=150, facecolor='#1a1a2e', edgecolor='none', bbox_inches='tight')
buf.seek(0)
img_base64 = base64.b64encode(buf.read()).decode('utf-8')
plt.close()

print("\n" + "="*70)
print("NUMERICAL ANALYSIS")
print("="*70)
print("\nEnergy in each snapshot (∫u²dx):")
for i, (t, _) in enumerate(snapshots):
    print(f"  t = {t:.3f}: E = {energies[i]:.6f}")
print(f"\nEnergy variation: {100*(max(energies)-min(energies))/energies[0]:.2f}%")
print("\n(Note: Energy not conserved due to reflecting boundary conditions)")

print("\n" + "="*70)
print("CONNECTION TO NUMERICAL RELATIVITY")
print("="*70)
print("\nKey concepts shared with full numerical relativity:")
print("\n1. BSSN evolution equations are also hyperbolic (wave-like)")
print("2. CFL condition ensures stability: $ \\Delta t < \\Delta x / c $")
print("3. Finite difference methods discretize spacetime")
print("4. Constraint preservation requires careful formulation")
print("\nIn full NR, we evolve:")
print("• Conformal metric: $ \\tilde{\\gamma}_{ij} $")
print("• Conformal factor: $ \\phi $")
print("• Traceless extrinsic curvature: $ \\tilde{A}_{ij} $")
print("• Mean curvature: $ K $")

print(f"\n<img src='data:image/png;base64,{img_base64}' alt='Wave Equation Evolution' style='max-width:100%;'/>")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran: BSSN Variable Computation

This program computes the key BSSN variables for Schwarzschild initial data, demonstrating the conformal decomposition and gauge conditions used in modern numerical relativity codes.

BSSN Variables Visualization

Python

Lapse function, conformal factor, and metric components vs radius

bssn_plot.py131 lines

import subprocess, tempfile, os, io, base64
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

fortran_source = r"""
program bssn_vars
    implicit none
    real(8), parameter :: G = 6.67430d-11
    real(8), parameter :: c = 2.99792458d8
    real(8), parameter :: M_sun = 1.98892d30
    real(8) :: M, r_s, r, alpha, psi4, chi_var, f
    integer :: i, n_pts

M = 10.0d0 * M_sun
    r_s = 2.0d0 * G * M / (c * c)

n_pts = 200
    print *, 'DATA_START'
    do i = 1, n_pts
        r = 2.01d0 * r_s + dble(i-1) * 48.0d0 * r_s / dble(n_pts-1)
        f = 1.0d0 - r_s / r
        alpha = sqrt(f)
        psi4 = (1.0d0 + r_s / (4.0d0 * r))**4
        chi_var = 1.0d0 / psi4
        write(*,'(F14.8,A,F14.8,A,F14.8,A,F14.8,A,F14.8)') &
            r/r_s, ',', alpha, ',', psi4, ',', chi_var, ',', 1.0d0/f
    end do
    print *, 'DATA_END'
end program bssn_vars
"""

workdir = tempfile.mkdtemp()
src = os.path.join(workdir, 'bssn.f90')
exe = os.path.join(workdir, 'bssn')

with open(src, 'w') as f:
    f.write(fortran_source)

comp = subprocess.run(['gfortran', '-O2', '-o', exe, src],
                      capture_output=True, text=True, timeout=30)
if comp.returncode != 0:
    print("Compilation error:", comp.stderr)
else:
    result = subprocess.run([exe], capture_output=True, text=True, timeout=30, cwd=workdir)
    output = result.stdout

lines = []
    in_data = False
    for line in output.split('\n'):
        if 'DATA_START' in line:
            in_data = True
            continue
        if 'DATA_END' in line:
            break
        if in_data and line.strip():
            lines.append(line.strip())

r_arr, alpha_arr, psi4_arr, chi_arr, grr_arr = [], [], [], [], []
    for line in lines:
        vals = [float(x) for x in line.split(',')]
        r_arr.append(vals[0])
        alpha_arr.append(vals[1])
        psi4_arr.append(vals[2])
        chi_arr.append(vals[3])
        grr_arr.append(vals[4])

r_arr = np.array(r_arr)
    alpha_arr = np.array(alpha_arr)
    psi4_arr = np.array(psi4_arr)
    chi_arr = np.array(chi_arr)
    grr_arr = np.array(grr_arr)

fig, axes = plt.subplots(1, 3, figsize=(18, 6))
    fig.patch.set_facecolor('#0f172a')
    fig.suptitle('BSSN Variables for Schwarzschild Initial Data', color='white', fontsize=14, fontweight='bold')

ax = axes[0]
    ax.set_facecolor('#1e293b')
    ax.plot(r_arr, alpha_arr, color='#22d3ee', linewidth=2, label='$\\alpha$ (lapse)')
    ax.axhline(y=1.0, color='#94a3b8', linestyle=':', alpha=0.5)
    ax.axvline(x=2.0, color='#ef4444', linestyle='--', alpha=0.7, label='Horizon ($r=r_s$)')
    ax.set_xlabel('$r / r_s$', color='#cbd5e1')
    ax.set_ylabel('Lapse function $\\alpha$', color='#cbd5e1')
    ax.set_title('Lapse: $\\alpha = \\sqrt{1 - r_s/r}$', color='white')
    ax.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')
    ax.tick_params(colors='#cbd5e1')
    for spine in ax.spines.values():
        spine.set_color('#334155')
    ax.grid(True, color='#334155', alpha=0.4)
    ax.set_ylim(0, 1.1)

ax2 = axes[1]
    ax2.set_facecolor('#1e293b')
    ax2.plot(r_arr, psi4_arr, color='#f472b6', linewidth=2, label='$\\psi^4$ (conformal)')
    ax2.plot(r_arr, chi_arr, color='#4ade80', linewidth=2, label='$\\chi = \\psi^{-4}$')
    ax2.axvline(x=2.0, color='#ef4444', linestyle='--', alpha=0.7)
    ax2.set_xlabel('$r / r_s$', color='#cbd5e1')
    ax2.set_ylabel('Conformal factor', color='#cbd5e1')
    ax2.set_title('Conformal Decomposition', color='white')
    ax2.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')
    ax2.tick_params(colors='#cbd5e1')
    for spine in ax2.spines.values():
        spine.set_color('#334155')
    ax2.grid(True, color='#334155', alpha=0.4)

ax3 = axes[2]
    ax3.set_facecolor('#1e293b')
    ax3.plot(r_arr, grr_arr, color='#fbbf24', linewidth=2, label='$g_{rr} = (1-r_s/r)^{-1}$')
    ax3.axvline(x=2.0, color='#ef4444', linestyle='--', alpha=0.7, label='Horizon')
    ax3.set_xlabel('$r / r_s$', color='#cbd5e1')
    ax3.set_ylabel('$g_{rr}$', color='#cbd5e1')
    ax3.set_title('Radial Metric Component', color='white')
    ax3.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1')
    ax3.tick_params(colors='#cbd5e1')
    for spine in ax3.spines.values():
        spine.set_color('#334155')
    ax3.grid(True, color='#334155', alpha=0.4)
    ax3.set_ylim(0, 20)

plt.tight_layout()
    buf = io.BytesIO()
    plt.savefig(buf, format='png', dpi=150, bbox_inches='tight', facecolor=fig.get_facecolor())
    buf.seek(0)
    img_b64 = base64.b64encode(buf.read()).decode()
    buf.close()
    plt.close()
    print("BSSN initial data: K=0 (time-symmetric), lapse -> 0 at horizon")
    print(f'<img src="data:image/png;base64,{img_b64}" />')

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← ADM Formalism Next: Gravitational Lensing →