on Billiards

Overview

This post examines the problem of speed reconstruction in a heterogeneous 2D disk based on Hamiltonian ray dynamics. We explore the linearization of the scattering relation, the resulting linear systems, and the “Foliation” algorithm for stable inward reconstruction.

🏷️ Hamiltonian Ray Dynamics

On a heterogeneous 2D disk, we consider geodesics in phase space governed by a Hamiltonian of the form $H(x,\xi )=c^2(x)|\xi {|}^2$, which describes the motion of a point billiard. To fully recover the speed profile $c(x)$, we require the scattering relation $(X_s,X_r,t_{sr})$, where $X_s$ is the input phase coordinate, $X_r$ is the output phase coordinate, and $t_{sr}$ is the total travel time.

In practice, we deploy equi-spaced sources on the boundary circle with equi-spaced input angles. The resulting data set is represented by a matrix where each row corresponds to a single ray’s trajectory and boundary measurements.

🏷️ Linearization and Reconstruction

The reconstruction of $c(x)$ via linearization leads to a system of time-integral equations. The ray trajectory $X(s)$ satisfies:

$$\frac{dX}{ds}=L(X),X(0)=X_s,X(t_{sr})=X_r$$

where $L=\mathrm{diag}(1,-1)\partial H/\partial X$. The Jacobian $J(s)=\partial X(s)/\partial X(0)$ evolves according to:

$$\frac{dJ}{ds}=M(X)J,J(0)=I$$

Applying the principle of variations, the mismatch in the boundary coordinates is related to the speed perturbation $\delta c$ by:

$$X_r-X_s\simeq {\int }_0^{t_{sr}}J(t_{sr})J(t{)}^{-1}S(\delta c,X(t))dt$$

Discretizing the integral using the trapezoid rule over a spatial grid yields a large, sparse linear system $D\simeq A\delta c$. In this formulation, the data vector $D$ represents the residuals or observed mismatch in output phase coordinates $X_r$ between the heterogeneous medium and the background reference. The sensitivity matrix $A$ is a sparse Jacobian where each row encodes the integrated influence of a single ray’s trajectory on local speed perturbations. Our objective is to solve for $\delta c$, the unknown vector of speed perturbations at each grid point that describes the deviation from the reference speed profile.

🏷️ Foliation

The matrix $A$ exhibits a rank ordering property: rays that stay near the boundary involve fewer unknown grid points than rays that penetrate the interior. By reordering the measurements according to the number of non-zero entries in each row (the ray’s footprint), we can solve the system iteratively from the boundary inward (Zhao & Zhong, 2019).

The Foliation Algorithm

1. Identify known boundary values and move their indices to a “solved” set $S_2$.
2. Rearrange the remaining equations in $A$ by their sparsity (nnz count).
3. Solve the subset of equations that depend primarily on $S_2$ and a small layer of unknown points.
4. Update the speed profile, update $S_2$, and rerun the forward problem to generate new residuals.
5. Repeat until the interior is fully resolved.

🏷️ Complexity and Stability

The linearized scheme provides at most first-order convergence. Computation is dominated by the numerical integration of the Hamiltonian system, which is 4-dimensional per ray. For large-scale measurements ($>10^3$), parallel solvers are essential.

The exponential accumulation of error (typical of time-integrated dynamics) makes fine-grid reconstruction numerically unstable without proper regularization or high-density boundary sampling.

📊 Numerical Verification

• Simulation Code: content/codes/2017 Fall/billiards_reconstruction.py

The foliation algorithm was verified through an iterative Gauss-Newton Confidence Propagation scheme.

Complex Medium Reconstruction

Fig 1: Left: True speed perturbation $\delta c$. Middle: Reconstructed $\delta c$ at 20th iteration. Right: Fidelity map showing trust distribution.

🏷️ Performance Analysis

• Measurement Density: 2400 rays from 60 boundary sensors.
• Ray Dynamics: Hamiltonian integration via solve_ivp with event-based boundary detection.
• Stability: Active pixel confidence reached $>0.9$, confirming the robustness of the fidelity-weighted update scheme.

🔗 See Also

• on Adjoint Framework --- Inverse problems often utilize adjoint state methods for gradient computation, providing a scalable alternative to the linearized Jacobian approach used here.
• on Faster Codes I --- The parallelization strategy used in the ray-tracing simulation reflects the performance optimization principles discussed for Python-based scientific computing.

📚 References

🐻  Zhao, H. & Zhong, Y. 2019. A hybrid adaptive phase space method for reflection traveltime tomography. SIAM Journal on Imaging Sciences 12(1), 28–53.