Benchmarks

ReferenceSolution dataclass

ReferenceSolution(x, t, u)

A solution sampled on a tensor-product (x, t) grid.

Attributes:

Name	Type	Description
x	ndarray	Spatial grid, shape (nx,).
t	ndarray	Temporal grid, shape (nt,).
u	ndarray	Solution values, shape (nx, nt) with u[i, j] = u(x[i], t[j]).

flat_points

flat_points()

Grid coordinates flattened to (nx * nt, 1) columns, ij order.

Source code in pinnrl/benchmarks/reference.py

def flat_points(self) -> tuple[np.ndarray, np.ndarray]:
    """Grid coordinates flattened to ``(nx * nt, 1)`` columns, ``ij`` order."""
    xx, tt = np.meshgrid(self.x, self.t, indexing="ij")
    return xx.reshape(-1, 1), tt.reshape(-1, 1)

relative_l2

relative_l2(u_pred)

||u_pred - u||_2 / ||u||_2 over the grid.

Source code in pinnrl/benchmarks/reference.py

def relative_l2(self, u_pred: np.ndarray) -> float:
    """``||u_pred - u||_2 / ||u||_2`` over the grid."""
    u_ref = self.u.reshape(-1)
    diff = np.asarray(u_pred).reshape(-1) - u_ref
    return float(np.linalg.norm(diff) / np.linalg.norm(u_ref))

burgers_reference

burgers_reference(
    nu=0.01 / np.pi,
    nx=256,
    nt=101,
    t_max=1.0,
    n_quad=6001,
    z_max=15.0,
)

Burgers u_t + u u_x = nu u_xx on [-1, 1] with u0 = -sin(pi x).

Cole-Hopf gives

u(x, t) = -∫ sin(pi (x - η)) f(x - η) G(η, t) dη / ∫ f(x - η) G(η, t) dη

with f(y) = exp(-cos(pi y) / (2 pi nu)) and the heat kernel G. Substituting η = sqrt(4 nu t) z turns G into exp(-z²); the integrand is periodic times Gaussian, so the trapezoid rule converges spectrally. z_max = 15 keeps the far peaks of f (which dominate near x = 0 once the shock forms) inside the window.

Parameters:

Name	Type	Description	Default
nu	float	Viscosity. The default is the classic 0.01 / pi benchmark.	0.01 / pi
nx	int	Number of spatial grid points on [-1, 1].	256
nt	int	Number of time grid points on [0, t_max].	101
t_max	float	Final time.	1.0
n_quad	int	Trapezoid nodes on [-z_max, z_max].	6001
z_max	float	Half-width of the quadrature window in z.	15.0

Returns:

Type	Description
ReferenceSolution	The reference field on the (nx, nt) grid.

Source code in pinnrl/benchmarks/reference.py

def burgers_reference(
    nu: float = 0.01 / np.pi,
    nx: int = 256,
    nt: int = 101,
    t_max: float = 1.0,
    n_quad: int = 6001,
    z_max: float = 15.0,
) -> ReferenceSolution:
    """Burgers ``u_t + u u_x = nu u_xx`` on ``[-1, 1]`` with ``u0 = -sin(pi x)``.

    Cole-Hopf gives

        u(x, t) = -∫ sin(pi (x - η)) f(x - η) G(η, t) dη / ∫ f(x - η) G(η, t) dη

    with ``f(y) = exp(-cos(pi y) / (2 pi nu))`` and the heat kernel ``G``.
    Substituting ``η = sqrt(4 nu t) z`` turns ``G`` into ``exp(-z²)``; the
    integrand is periodic times Gaussian, so the trapezoid rule converges
    spectrally. ``z_max = 15`` keeps the far peaks of ``f`` (which dominate
    near ``x = 0`` once the shock forms) inside the window.

    Args:
        nu: Viscosity. The default is the classic ``0.01 / pi`` benchmark.
        nx: Number of spatial grid points on ``[-1, 1]``.
        nt: Number of time grid points on ``[0, t_max]``.
        t_max: Final time.
        n_quad: Trapezoid nodes on ``[-z_max, z_max]``.
        z_max: Half-width of the quadrature window in ``z``.

    Returns:
        The reference field on the ``(nx, nt)`` grid.
    """
    x = np.linspace(-1.0, 1.0, nx)
    t = np.linspace(0.0, t_max, nt)
    u = burgers_cole_hopf(x[:, None], t[None, :], nu, n_quad=n_quad, z_max=z_max)
    return ReferenceSolution(x=x, t=t, u=u)

allen_cahn_reference

allen_cahn_reference(
    epsilon=np.sqrt(0.0002),
    t_max=5.0,
    nx=257,
    nt=101,
    refine=8,
    boundary_value=-1.0,
)

Allen-Cahn u_t = eps² u_xx + u - u³ on [-1, 1], u0 = x² cos(pi x).

Dirichlet u(±1, t) = boundary_value. The defaults are the u_t = 0.001 u_xx + 5 (u - u³) benchmark on t ∈ [0, 1] rewritten in the time variable τ = 5 t, which matches the reaction coefficient of :class:pinnrl.pdes.allen_cahn.AllenCahnEquation.

The solver grid is refine times finer than the output grid and the output points are a subset of it, so no interpolation error enters.

Parameters:

Name	Type	Description	Default
epsilon	float	Interface parameter in the residual u_t - eps² u_xx - u + u³.	sqrt(0.0002)
t_max	float	Final time.	5.0
nx	int	Output spatial points on [-1, 1].	257
nt	int	Output time points on [0, t_max].	101
refine	int	Solver-grid refinement factor relative to the output grid.	8
boundary_value	float	Dirichlet value at both ends.	-1.0

Returns:

Type	Description
ReferenceSolution	The reference field on the (nx, nt) grid.

Source code in pinnrl/benchmarks/reference.py

def allen_cahn_reference(
    epsilon: float = np.sqrt(0.0002),
    t_max: float = 5.0,
    nx: int = 257,
    nt: int = 101,
    refine: int = 8,
    boundary_value: float = -1.0,
) -> ReferenceSolution:
    """Allen-Cahn ``u_t = eps² u_xx + u - u³`` on ``[-1, 1]``, ``u0 = x² cos(pi x)``.

    Dirichlet ``u(±1, t) = boundary_value``. The defaults are the
    ``u_t = 0.001 u_xx + 5 (u - u³)`` benchmark on ``t ∈ [0, 1]`` rewritten
    in the time variable ``τ = 5 t``, which matches the reaction coefficient
    of :class:`pinnrl.pdes.allen_cahn.AllenCahnEquation`.

    The solver grid is ``refine`` times finer than the output grid and the
    output points are a subset of it, so no interpolation error enters.

    Args:
        epsilon: Interface parameter in the residual ``u_t - eps² u_xx - u + u³``.
        t_max: Final time.
        nx: Output spatial points on ``[-1, 1]``.
        nt: Output time points on ``[0, t_max]``.
        refine: Solver-grid refinement factor relative to the output grid.
        boundary_value: Dirichlet value at both ends.

    Returns:
        The reference field on the ``(nx, nt)`` grid.
    """
    n_fine = (nx - 1) * refine + 1
    x_fine = np.linspace(-1.0, 1.0, n_fine)
    dx = x_fine[1] - x_fine[0]
    t = np.linspace(0.0, t_max, nt)
    eps2 = epsilon**2

    u0 = x_fine**2 * np.cos(np.pi * x_fine)

    def rhs(_t: float, v: np.ndarray) -> np.ndarray:
        full = np.concatenate(([boundary_value], v, [boundary_value]))
        u_xx = (full[2:] - 2.0 * full[1:-1] + full[:-2]) / dx**2
        return eps2 * u_xx + v - v**3

    sol = solve_ivp(rhs, (0.0, t_max), u0[1:-1], t_eval=t, method="RK45", rtol=1e-9, atol=1e-11)
    if not sol.success:
        raise RuntimeError(f"Allen-Cahn reference solve failed: {sol.message}")

    u_fine = np.empty((n_fine, nt))
    u_fine[0, :] = boundary_value
    u_fine[-1, :] = boundary_value
    u_fine[1:-1, :] = sol.y
    return ReferenceSolution(x=x_fine[::refine], t=t, u=u_fine[::refine, :])

StudySettings dataclass

StudySettings(
    num_steps=15000,
    n_collocation=1000,
    resample_every=100,
    learning_rate=0.001,
    hidden_dims=(64, 64, 64, 64),
    n_checkpoints=50,
    rad_k=1.0,
    rad_c=1.0,
    rad_pool_factor=20,
)

Shared knobs for every run in a study.

Attributes:

Name	Type	Description
num_steps	int	Optimizer steps per run (the equal training budget).
n_collocation	int	Collocation points per step.
resample_every	int	Steps between resampling for uniform_periodic, rad and rl_bandit.
learning_rate	float	Adam learning rate (constant).
hidden_dims	tuple[int, ...]	Widths of the tanh MLP.
n_checkpoints	int	Evaluations along the run; num_steps / n_checkpoints must be a multiple of resample_every.
rad_k	float	RAD exponent k.
rad_c	float	RAD uniform floor c.
rad_pool_factor	int	RAD candidate pool size as a multiple of n_collocation.

chunk property

chunk

Steps between checkpoints.

run_sampling_study

run_sampling_study(
    pdes=STUDY_PDES,
    strategies=STUDY_STRATEGIES,
    seeds=tuple(range(10)),
    settings=None,
    workers=1,
    on_result=None,
)

Run every (pde, strategy, seed) combination.

Tasks are ordered by PDE, then seed, then strategy, so the strategies for one seed run side by side under similar machine load. workers > 1 uses a spawn-based process pool; workers = 1 runs inline.

Parameters:

Name	Type	Description	Default
pdes	Sequence[str]	Subset of :data:STUDY_PDES.	STUDY_PDES
strategies	Sequence[str]	Subset of :data:STUDY_STRATEGIES.	STUDY_STRATEGIES
seeds	Sequence[int]	Seeds to run; each fixes the initial weights.	tuple(range(10))
settings	StudySettings | None	Shared study settings.	None
workers	int	Parallel worker processes.	1
on_result	Callable[[StudyRun], None] | None	Optional callback invoked as each run finishes.	None

Returns:

Name	Type	Description
One	list[StudyRun]	class:StudyRun per combination, in task order.

Source code in pinnrl/benchmarks/sampling_study.py

def run_sampling_study(
    pdes: Sequence[str] = STUDY_PDES,
    strategies: Sequence[str] = STUDY_STRATEGIES,
    seeds: Sequence[int] = tuple(range(10)),
    settings: StudySettings | None = None,
    workers: int = 1,
    on_result: Callable[[StudyRun], None] | None = None,
) -> list[StudyRun]:
    """Run every ``(pde, strategy, seed)`` combination.

    Tasks are ordered by PDE, then seed, then strategy, so the strategies for
    one seed run side by side under similar machine load. ``workers > 1`` uses a spawn-based
    process pool; ``workers = 1`` runs inline.

    Args:
        pdes: Subset of :data:`STUDY_PDES`.
        strategies: Subset of :data:`STUDY_STRATEGIES`.
        seeds: Seeds to run; each fixes the initial weights.
        settings: Shared study settings.
        workers: Parallel worker processes.
        on_result: Optional callback invoked as each run finishes.

    Returns:
        One :class:`StudyRun` per combination, in task order.
    """
    settings = settings or StudySettings()
    for name in pdes:
        if name not in _PDE_BUILDERS:
            raise ValueError(f"Unknown pde {name!r}. Supported: {STUDY_PDES}")
    for name in strategies:
        if name not in STUDY_STRATEGIES:
            raise ValueError(f"Unknown strategy {name!r}. Supported: {STUDY_STRATEGIES}")

    references = {name: build_reference(name) for name in pdes}
    tasks = [
        (pde, strategy, seed, settings, references[pde])
        for pde in pdes
        for seed in seeds
        for strategy in strategies
    ]

    results: list[StudyRun] = []
    if workers <= 1:
        for task in tasks:
            run = run_study_case(*task)
            results.append(run)
            if on_result:
                on_result(run)
        return results

    with ProcessPoolExecutor(max_workers=workers, mp_context=mp.get_context("spawn")) as pool:
        for run in pool.map(_worker, tasks):
            results.append(run)
            if on_result:
                on_result(run)
    return results

run_study_case

run_study_case(
    pde_name,
    strategy,
    seed,
    settings,
    reference,
    device=None,
)

Train one PINN for settings.num_steps steps under one sampler.

Source code in pinnrl/benchmarks/sampling_study.py

def run_study_case(
    pde_name: str,
    strategy: str,
    seed: int,
    settings: StudySettings,
    reference: ReferenceSolution,
    device: torch.device | None = None,
) -> StudyRun:
    """Train one PINN for ``settings.num_steps`` steps under one sampler."""
    if pde_name not in _PDE_BUILDERS:
        raise ValueError(f"Unknown pde {pde_name!r}. Supported: {STUDY_PDES}")
    if strategy not in STUDY_STRATEGIES:
        raise ValueError(f"Unknown strategy {strategy!r}. Supported: {STUDY_STRATEGIES}")
    chunk = settings.chunk
    if chunk * settings.n_checkpoints != settings.num_steps or chunk % settings.resample_every:
        raise ValueError(
            "num_steps must split into n_checkpoints chunks that are each a multiple "
            f"of resample_every (got {settings.num_steps}, {settings.n_checkpoints}, "
            f"{settings.resample_every})."
        )
    device = device or torch.device("cpu")

    torch.manual_seed(seed)
    np.random.seed(seed)
    pde = _PDE_BUILDERS[pde_name](device)
    trainer = _build_trainer(pde, settings, device)
    torch.manual_seed(seed + 10_000)  # sampling stream, independent of the init

    overhead = {"evals": 0}
    progress = {"value": 0.0}
    sampler = _make_sampler(strategy, pde, settings, device, overhead)
    resample_every = 1 if strategy in _EVERY_STEP else settings.resample_every

    def sample_fn(model, _chunk_progress):
        return sampler(model, progress["value"])

    run = StudyRun(pde_name, strategy, seed, float("nan"), 0, 0.0, 0)
    for _ in range(settings.n_checkpoints):
        out = trainer.train_steps(
            num_steps=chunk, sample_fn=sample_fn, resample_every=resample_every
        )
        run.steps += out["steps"]
        run.train_time_s += out["train_time_s"]
        progress["value"] = run.steps / settings.num_steps
        run.residual_evals = run.steps * settings.n_collocation + overhead["evals"]
        run.curve_step.append(run.steps)
        run.curve_time_s.append(run.train_time_s)
        run.curve_evals.append(run.residual_evals)
        run.curve_rel_l2.append(_relative_l2(trainer.model, reference, device))
    run.rel_l2 = run.curve_rel_l2[-1]
    return run

summarize

summarize(runs, budget='steps')

Median/IQR per (pde, strategy) plus a paired comparison against uniform.

budget="steps" compares final errors (equal training points). budget="evals" compares each run at the last checkpoint whose residual evaluations (training plus sampler overhead) do not exceed the smallest total spent by any run on that PDE.

Source code in pinnrl/benchmarks/sampling_study.py

def summarize(runs: Sequence[StudyRun], budget: str = "steps") -> list[StudySummary]:
    """Median/IQR per ``(pde, strategy)`` plus a paired comparison against ``uniform``.

    ``budget="steps"`` compares final errors (equal training points).
    ``budget="evals"`` compares each run at the last checkpoint whose residual
    evaluations (training plus sampler overhead) do not exceed the smallest
    total spent by any run on that PDE.
    """
    from scipy.stats import wilcoxon

    if budget not in ("steps", "evals"):
        raise ValueError(f"budget must be 'steps' or 'evals', got {budget!r}")
    cap: dict[str, int] = {}
    if budget == "evals":
        for run in runs:
            cap[run.pde] = min(cap.get(run.pde, run.residual_evals), run.residual_evals)

    by_key: dict[tuple[str, str], dict[int, tuple[float, float]]] = {}
    for run in runs:
        point = _at_budget(run, cap.get(run.pde))
        by_key.setdefault((run.pde, run.strategy), {})[run.seed] = point

    summaries: list[StudySummary] = []
    for (pde, strategy), seed_vals in sorted(by_key.items()):
        errors = np.array([v[0] for v in seed_vals.values()])
        baseline = by_key.get((pde, "uniform"), {})
        paired = sorted(set(seed_vals) & set(baseline))
        wins = sum(seed_vals[s][0] < baseline[s][0] for s in paired)
        p_value = float("nan")
        if strategy != "uniform" and len(paired) >= 5:
            diffs = [seed_vals[s][0] - baseline[s][0] for s in paired]
            if any(d != 0 for d in diffs):
                p_value = float(wilcoxon(diffs).pvalue)
        summaries.append(
            StudySummary(
                budget=budget,
                pde=pde,
                strategy=strategy,
                n=len(errors),
                median_rel_l2=float(np.median(errors)),
                q25_rel_l2=float(np.percentile(errors, 25)),
                q75_rel_l2=float(np.percentile(errors, 75)),
                median_time_s=float(np.median([v[1] for v in seed_vals.values()])),
                wins_vs_uniform=wins if strategy != "uniform" else 0,
                p_vs_uniform=p_value,
            )
        )
    return summaries

solve_heat_1d

solve_heat_1d(
    alpha=0.1,
    frequency=1.0,
    domain=(0.0, 1.0),
    t_max=1.0,
    nx=101,
    nt=4001,
)

Explicit-Euler FDM solver for the 1-D heat equation on a periodic domain.

Solves u_t = α·u_xx with initial condition sin(2π·f·x) and periodic boundaries. The defaults give a diffusion number r = α·dt/dx² ≈ 0.26, about half the explicit-stability cap of 0.5.

Parameters:

Name	Type	Description	Default
alpha	float	Thermal diffusivity.	0.1
frequency	float	Spatial frequency of the sinusoidal initial condition.	1.0
domain	tuple[float, float]	Spatial domain (x_min, x_max).	(0.0, 1.0)
t_max	float	Final time.	1.0
nx	int	Number of distinct points on the periodic grid (dx = (x_max - x_min) / nx).	101
nt	int	Number of time levels including t = 0 (nt - 1 steps).	4001

Returns:

Type	Description
FDMResult	class:FDMResult containing the full solution grid, the analytical
FDMResult	solution at t_max, and L2/max error metrics.

Raises:

Type	Description
ValueError	If the explicit-Euler scheme would be unstable for the chosen (nx, nt, alpha).

Source code in pinnrl/benchmarks/fdm.py

def solve_heat_1d(
    alpha: float = 0.1,
    frequency: float = 1.0,
    domain: tuple[float, float] = (0.0, 1.0),
    t_max: float = 1.0,
    nx: int = 101,
    nt: int = 4001,
) -> FDMResult:
    """Explicit-Euler FDM solver for the 1-D heat equation on a periodic domain.

    Solves ``u_t = α·u_xx`` with initial condition ``sin(2π·f·x)`` and
    periodic boundaries. The defaults give a diffusion number
    ``r = α·dt/dx² ≈ 0.26``, about half the explicit-stability cap of 0.5.

    Args:
        alpha: Thermal diffusivity.
        frequency: Spatial frequency of the sinusoidal initial condition.
        domain: Spatial domain ``(x_min, x_max)``.
        t_max: Final time.
        nx: Number of distinct points on the periodic grid (``dx = (x_max - x_min) / nx``).
        nt: Number of time levels including ``t = 0`` (``nt - 1`` steps).

    Returns:
        :class:`FDMResult` containing the full solution grid, the analytical
        solution at ``t_max``, and L2/max error metrics.

    Raises:
        ValueError: If the explicit-Euler scheme would be unstable for the
            chosen ``(nx, nt, alpha)``.
    """
    x_min, x_max = float(domain[0]), float(domain[1])
    # Periodic grid: x_max is the same point as x_min, so it is excluded.
    x = np.linspace(x_min, x_max, nx, endpoint=False, dtype=np.float64)
    t = np.linspace(0.0, t_max, nt, dtype=np.float64)
    dx = (x_max - x_min) / nx
    dt = t_max / (nt - 1)
    r = alpha * dt / (dx**2)
    if r > 0.5:
        raise ValueError(
            f"FDM heat scheme unstable: r = α·dt/dx² = {r:.3f} > 0.5. "
            f"Increase nt (currently {nt}) or decrease nx (currently {nx})."
        )

    u = np.zeros((nt, nx), dtype=np.float64)
    u[0] = np.sin(2.0 * np.pi * frequency * x)

    start = time.perf_counter()
    for n in range(nt - 1):
        # Periodic stencil via np.roll: u_xx ≈ (u[i+1] - 2u[i] + u[i-1]) / dx²
        lap = np.roll(u[n], -1) - 2.0 * u[n] + np.roll(u[n], 1)
        u[n + 1] = u[n] + r * lap
    wall = time.perf_counter() - start

    u_exact_final = _heat_exact(x, t_max, alpha, frequency)
    err = u[-1] - u_exact_final
    return FDMResult(
        u=u,
        u_exact_final=u_exact_final,
        x=x,
        t=t,
        l2_error=float(np.linalg.norm(err) / np.linalg.norm(u_exact_final)),
        mse=float(np.mean(err**2)),
        max_error=float(np.max(np.abs(err))),
        wall_time_s=wall,
    )

solve_wave_1d

solve_wave_1d(
    c=1.0, domain=(0.0, 1.0), t_max=1.0, nx=201, nt=4001
)

Leap-frog FDM solver for the 1-D wave equation on a periodic domain.

Solves u_tt = c²·u_xx with the travelling-wave initial condition u(x, 0) = sin(2πx), u_t(x, 0) = -2π·c·cos(2πx) so that the closed-form solution is the right-moving wave sin(2π(x − c·t)). This is not :class:pinnrl.pdes.wave_equation.WaveEquation's standing-wave exact solution, which starts at rest.

Parameters:

Name	Type	Description	Default
c	float	Wave speed.	1.0
domain	tuple[float, float]	Spatial domain (x_min, x_max).	(0.0, 1.0)
t_max	float	Final time.	1.0
nx	int	Number of distinct points on the periodic grid (dx = (x_max - x_min) / nx).	201
nt	int	Number of time levels including t = 0 (nt - 1 steps).	4001

Returns:

Type	Description
FDMResult	class:FDMResult containing the full solution grid, the analytical
FDMResult	solution at t_max, and L2/max error metrics.

Raises:

Type	Description
ValueError	If the leap-frog CFL condition c·dt/dx ≤ 1 is violated.

Source code in pinnrl/benchmarks/fdm.py

def solve_wave_1d(
    c: float = 1.0,
    domain: tuple[float, float] = (0.0, 1.0),
    t_max: float = 1.0,
    nx: int = 201,
    nt: int = 4001,
) -> FDMResult:
    """Leap-frog FDM solver for the 1-D wave equation on a periodic domain.

    Solves ``u_tt = c²·u_xx`` with the travelling-wave initial condition
    ``u(x, 0) = sin(2πx)``, ``u_t(x, 0) = -2π·c·cos(2πx)`` so that the
    closed-form solution is the right-moving wave ``sin(2π(x − c·t))``. This
    is not :class:`pinnrl.pdes.wave_equation.WaveEquation`'s standing-wave
    exact solution, which starts at rest.

    Args:
        c: Wave speed.
        domain: Spatial domain ``(x_min, x_max)``.
        t_max: Final time.
        nx: Number of distinct points on the periodic grid (``dx = (x_max - x_min) / nx``).
        nt: Number of time levels including ``t = 0`` (``nt - 1`` steps).

    Returns:
        :class:`FDMResult` containing the full solution grid, the analytical
        solution at ``t_max``, and L2/max error metrics.

    Raises:
        ValueError: If the leap-frog CFL condition ``c·dt/dx ≤ 1`` is violated.
    """
    x_min, x_max = float(domain[0]), float(domain[1])
    # Periodic grid: x_max is the same point as x_min, so it is excluded.
    x = np.linspace(x_min, x_max, nx, endpoint=False, dtype=np.float64)
    t = np.linspace(0.0, t_max, nt, dtype=np.float64)
    dx = (x_max - x_min) / nx
    dt = t_max / (nt - 1)
    cfl = c * dt / dx
    if cfl > 1.0:
        raise ValueError(
            f"FDM wave scheme violates CFL: c·dt/dx = {cfl:.3f} > 1. "
            f"Increase nt (currently {nt}) or decrease nx (currently {nx})."
        )

    u = np.zeros((nt, nx), dtype=np.float64)
    u[0] = np.sin(2.0 * np.pi * x)
    # First step from u_t(x, 0) = -2πc·cos(2πx). Taylor expansion:
    #   u(x, dt) ≈ u(x, 0) + dt·u_t + 0.5·dt²·u_tt   with u_tt = c²·u_xx
    u_t0 = -2.0 * np.pi * c * np.cos(2.0 * np.pi * x)
    lap0 = np.roll(u[0], -1) - 2.0 * u[0] + np.roll(u[0], 1)
    u[1] = u[0] + dt * u_t0 + 0.5 * (cfl**2) * lap0

    start = time.perf_counter()
    for n in range(1, nt - 1):
        lap = np.roll(u[n], -1) - 2.0 * u[n] + np.roll(u[n], 1)
        u[n + 1] = 2.0 * u[n] - u[n - 1] + (cfl**2) * lap
    wall = time.perf_counter() - start

    u_exact_final = _wave_exact(x, t_max, c)
    err = u[-1] - u_exact_final
    return FDMResult(
        u=u,
        u_exact_final=u_exact_final,
        x=x,
        t=t,
        l2_error=float(np.linalg.norm(err) / np.linalg.norm(u_exact_final)),
        mse=float(np.mean(err**2)),
        max_error=float(np.max(np.abs(err))),
        wall_time_s=wall,
    )