Long-Axis Split Method in PDEs

• Long-Axis Split Method is a numerical technique that decouples singular contributions along a principal axis in PDEs to facilitate standard discretization.
• It effectively addresses mixed-dimensional challenges in elliptic and wave propagation models, yielding optimal finite element convergence rates.
• The method leverages analytical singularity subtraction to overcome mesh bottlenecks in applications like vascular flow simulations and random media.

The long-axis split method, also called singularity subtraction or phase-screen Strang splitting (depending on context), refers to a class of splitting methods specialized for PDEs with strong anisotropy along a main axis. These include mixed-dimensional elliptic equations with 1D line sources embedded in 3D domains and paraxial or Itô–Schrödinger wave equations with principal propagation direction. The technique isolates analytically explicit, low-regularity contributions linked to the long axis, allowing numerical approximation of a regular correction via standard discretizations, thus overcoming severe singularity and mesh-resolution bottlenecks. The method is prominent in simulation of vascular flows and wave propagation in random media, enabling efficient and optimal-accuracy computations even for highly singular or random-inhomogeneous problems (Gjerde et al., 2018, Bal et al., 1 Mar 2025).

1. Formulation for Elliptic Equations with Line Sources

The long-axis split method was introduced for elliptic problems with Dirac measure sources supported on a network of line segments $\Lambda=\cup_{i=1}^n \Lambda_i$ in a domain $\Omega\subset\mathbb{R}^3$. The model is: $$-\nabla\cdot(\kappa\nabla u) = f \delta_\Lambda \quad \text{in }\Omega, \qquad u=u_D\quad\text{on }\partial\Omega,$$ with $\kappa\in W^{2,\infty}(\Omega)$ strictly positive. The singularity due to the 1D measure in 3D causes the solution $u$ to fail to belong to $H^1(\Omega)$ globally, but $u$ exhibits piecewise $H^2$ regularity in directions parallel to the line segments.

The splitting theorem [(Gjerde et al., 2018), Thm 3.1] yields

$$u(x) = \frac{1}{4\pi} \left( \sum_{i=1}^n \frac{E_i[f] G_i(x)}{\kappa(x)} + w(x) \right),$$

where $G_i(x)$ is the logarithmic potential along $\Lambda_i$, computed as

$$G_i(x) = \int_0^{L_i} \frac{dt}{\| x - (a_i + \tau_i t) \| } = \ln\left( \frac{ r_{b,i} + L_i + \tau_i \cdot (a_i - x) }{ r_{a,i} + \tau_i \cdot (a_i - x) } \right),$$

with $r_{a,i} = |x - a_i|$, $r_{b,i} = |x - b_i|$, and $E_i[f]$ is any smooth extension of $f_i$ off the line. The explicit sum $u_s$ encapsulates the singular field.

The regular correction $w\in H^1(\Omega)$ solves the standard elliptic PDE,

$$-\nabla\cdot (\kappa \nabla w) = F \quad\text{in }\Omega, \qquad w = w_D \text{ on } \partial\Omega,$$

where $F$ and $w_D$ are computed analytically from $E_i[f]$ and $G_i$.

2. Regularity and Anisotropy of the Singular Field

The explicit singular field $u_s$ fails to belong to $H^1(\Omega)$ due to its strong singularity near $\Lambda$, but is smooth and piecewise $H^2$ away from small tubular neighborhoods. Along the tangential direction of each segment, $u$ possesses the regularity of the source data $f$, i.e., $u_{\tau\tau}\in L^2(\Omega)$. This anisotropy is crucial for understanding error estimates and mesh requirements. The remainder $w$ inherits high regularity provided $E_i[f]$ is chosen such that $F\in L^2(\Omega)$, e.g., $E_i[f]$ is constant perpendicular to $\Lambda_i$ in a local tube.

Regularity summary:

Field	Global Sobolev Space	Tangential Direction Regularity
$u$	$W^{1,p}(\Omega)\ (p<2)$	$u_{\tau\tau}\in L^2(\Omega)$
$u\restriction_{\Omega\setminus U_\epsilon}$	$H^2$	$-$
$w$	$H^2(\Omega)$ (locally)	$-$

3. Numerical Approximation via Singularity Subtraction

The Galerkin finite element method (FEM) with singularity subtraction is implemented as:

1. Approximate $u_s(x)$ by interpolation onto the FE space (typically $P_1$ or $P_2$).
2. Assemble $F(x) = -\nabla\cdot(\kappa\nabla u_s^h)$ or use the closed-form analytic $F$.
3. Solve for $w_h$ in the FE space using the standard variational formulation.
4. Recover $u_h = u_s^h + w_h$.

Uniform tetrahedral/prism meshes are sufficient—no mesh refinement along line sources is required, since all singular effects are contained analytically in $u_s$. The global solve remains that of a standard 3D elliptic system with $L^2$ right-hand side.

4. Convergence and Computational Properties

Direct FEM on $u$ (discretizing the line source directly) converges sub-optimally and does not yield $H^1$-convergence due to global lack of $H^1$-regularity. The long-axis split/SSB-FEM approach achieves optimal rates on uniform meshes:

• $\|u - u_h\|_{L^2} = O(h^2)$,
• $\|u - u_h\|_{H^1} = O(h^1)$,

even when $u\notin H^1(\Omega)$, as $w$ is smooth and $u_s^h$ is known analytically.

For large-scale problems with $n\sim 3000$ line segments, e.g., brain vascular datasets, cut-cell integration is not necessary and all singular contributions reduce to closed-form $G_i(x)$. Efficient parallelization is feasible by vectorizing the evaluation of $G_i(x)E_i[f]$ and employing fast multipole methods as needed. The computational cost scales as $O(N_{\text{nodes}} + n \cdot N_q)$, where $N_q$ is the number of quadrature points (Gjerde et al., 2018).

5. Application to Paraxial and Stochastic Wave Propagation

In paraxial and Itô–Schrödinger wave propagation models, the long-axis split—termed phase-screen Strang splitting—exploits the natural separation between the main (long) axis ($z$) and transverse variables ($x$). The equations are:

• Paraxial: $$\partial_z u^\theta = i\kappa_1(z)\Delta_x u^\theta + i\kappa_2(z)\nu^\theta(z,x)u^\theta$$,
• Itô–Schrödinger: $$du = i\kappa_1(z)\Delta_x u\,dz - \frac12 \kappa_2(z)^2 R(0)u\,dz + i\kappa_2(z)u\,dB(z,x)$$,

with random or stochastic terms in the main axis (Bal et al., 1 Mar 2025).

The centered Strang-splitting scheme alternates half-step “phase-screen” potential propagations (in physical space) with full-step free propagation (in Fourier space) along $z$. The transverse Laplacian is diagonalized spectrally. Theoretical analysis shows mean-square pathwise convergence of order $\Delta z$ and moment error of order $(\Delta z)^2$ for the Strang case, even when the transverse (wavelength) scale is much smaller than the main axis discretization.

6. Algorithmic Realization

In both elliptic and wave propagation settings, the numerical method decouples the singular behavior from the regular correction, allowing the use of standard discretizations throughout. For the paraxial/Itô–Schrödinger model, the main algorithmic steps are:

• Initialize the field on the spatial grid.
• Precompute Laplacian symbols and Fourier grid for efficient spectral propagation.
• For each main-axis step:
  • Apply a half-screen (potential) step in physical space.
  • Evolve via free propagation using FFT/IFFT in Fourier space.
  • Apply the remaining half-screen step.
• In stochastic settings, replace deterministic steps with sampled stochastic increments as prescribed.

No mesh adaptation in the long-axis direction is required beyond uniform collocation, and all singular/random contributions are handled analytically or via explicit quadrature.

7. Impact, Advantages, and Large-Scale Applications

The long-axis split method enables efficient, accurate, and parallelizable simulation of mixed-dimensional and strongly anisotropic problems previously hindered by singularities and mesh constraints. Principal advantages include:

• Reduction of mixed-dimensional PDEs to standard elliptic/wave equations for smooth corrections.
• Elimination of mesh grading or cut-cell integration near singular sources.
• Optimal convergence and error rates on uniform meshes.
• Suitability for embarrassingly parallel computation of thousands of sources/interactions (Gjerde et al., 2018, Bal et al., 1 Mar 2025).

Demonstrations on brain-vascular networks ($\sim$3,000 segments) and wave propagation in random media confirm robust applicability, with runtime not degraded by singularity geometry or high source count. The method has established itself as a standard approach in challenging mixed-dimensional PDE settings requiring both mathematical rigor and computational scalability.