---
title: Long-Axis Split Method in PDEs
url: https://www.emergentmind.com/topics/long-axis-split-method
type: topic
---

# Long-Axis Split Method in PDEs

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 [1810.12979, 2503.00633].

## 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 [1810.12979, 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 [1810.12979].

## 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 [2503.00633].

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 [1810.12979, 2503.00633].

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.

Source: https://www.emergentmind.com/topics/long-axis-split-method