stencil_order.md

Stencil Order Verification

Claim

• The spatial operator in ref/code.cpp uses a 4th‑order central approximation of the second derivative per axis, with coefficients: [-1/12, 4/3, -2.5, 4/3, -1/12] at offsets {−2, −1, 0, +1, +2}.

Evidence

1. Coefficient match: the code applies

   • (-8.33333333e-2) = −1/12 to (x±2), (y±2), (z±2)
   • 1.333333330 = 4/3 to (x±1), (y±1), (z±1)
   • central −2.5 combined via r5 This is the standard 4th‑order stencil for ∂²/∂x² etc., then scaled by r2,r3,r4 = 1/h².

2. Convergence test:

   • analysis/verify_order.py computes the 3D Laplacian of f(x,y,z)=sin(axx)+sin(byy)+sin(czz) using this stencil and compares to the analytic laplacian.
   • It sweeps grid spacing h and fits the slope of log(error) vs log(h); result should be ≈ 4.

Run (two paths)

• Analytic reference (Python): python analysis/verify_order.py → analysis/stencil_convergence.png (slope ~4 in double‑precision arithmetic).
• Original C (FP32/FP64):
  • cd cpu_bench && make verify_order && python order_sweep.py --exe ./verify_order
  • make verify_order_fp64 && python order_sweep.py --exe ./verify_order_fp64
  • Plots: cpu_bench/stencil_convergence_sweep.png, cpu_bench/stencil_convergence_cpp_fp64.png.
  • Note: single‑step identity and finite precision yield an error floor ∝ ε/h², producing apparent slopes near −2 on practical meshes. Lower frequencies/coarser h push curves toward the 4th‑order asymptote.