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.

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